Sex Differences

Code

# load packages 
library(tidyverse)
library(gtsummary)
library(gt)
library(tidymodels)
library(censored)
library(paletteer)
library(survey)
library(patchwork)
theme_set(theme_light())
# paletteer_d("colorBlindness::PairedColor12Steps")
col1 = "#FF9932"; col2 = "#65CCFF"
col1 = "#CC5151"; col2 = "#422CB2"
# paletteer_d("ggthemes::colorblind")
col_vec = c("#000000FF", "#009E73FF", "#CC79A7FF", "#E69F00FF", "#D55E00FF", "#56B4E9FF", "#0072B2FF")

pa_old = readRDS(here::here("data", "accelerometry", 
                              "summarized", "pa_summary_AC_C_D.rds")) %>% 
  mutate(SEQN = as.character(SEQN)) %>% 
  group_by(SEQN) %>%
  summarize(across(-PAXDAY, ~mean(.x, na.rm = TRUE))) 

days = readRDS(here::here("data", "accelerometry", "inclusion_c_d.rds"))

wear = 
  days %>% 
  mutate(non_zero_vals = non_na_vals - zero_vals) %>% 
  filter(include_day) %>% 
  group_by(SEQN) %>% 
  summarize(wear_mins = mean(non_zero_vals)) %>% 
  mutate(SEQN = as.character(SEQN))

demo = readRDS(here::here("data", "demographics", "processed", "subset_C_D_tidy.rds")) 

NHANES 2003-2006

Survey Weighted Physical Activity Summaries by Sex

Code

df_small = pa_old %>% 
  left_join(demo) %>% 
  left_join(wear) %>% 
  filter(age_in_years_at_screening >= 18) %>%
  select(SEQN,
         data_release_cycle,
         gender,
         wear_mins,
         masked_variance_pseudo_psu,
         masked_variance_pseudo_stratum,
         full_sample_2_year_interview_weight,
         full_sample_2_year_mec_exam_weight,
         age_in_years_at_screening,
         total, 
         mean,
         mvpa, 
         mod,
         vig,
         mvpa_bout,
         vig_bout,
         starts_with("x"),
         missing_ac = na,
         median,
         logmean, 
         logmedian, 
         logtotal) 
labs = c("SEQN", "Data release cycle", "Gender",  "Wear minutes",
         "Pseudo PSU",
         "Psueudo stratum",
         "2 yr int weight", "2 yr exam weight",
         "Age", "Total AC", "Mean AC", "MVPA minutes",  "Moderate minutes", "Vigorous minutes",
         "MVPA (10 min bout)",
         "Vig (10 min bout)",
         "Time <5", "Time [5,10]","Time (10,25]","Time (25,50]", "Time (50,75]", "Time (75,100]", "Time (100,150]", "Time (150,200]", "Time (200,500]", 
         "Time (500,1000]", "Time (1000,2000]", "Time (2000,4000]","Time (4000,6000]","Time >6000", "Missing AC",
         "Median AC", "log10 mean AC", "log10 median AC", "log10 total AC") 

names(labs) = colnames(df_small)

df_small =
  df_small %>%
  labelled::set_variable_labels(!!!labs)

df_svy =
  df_small %>% 
  mutate(WTMEC4YR = full_sample_2_year_mec_exam_weight/2,
         WTMEC4YR_norm = WTMEC4YR/mean(WTMEC4YR, na.rm = TRUE)) %>%
  svydesign(ids = ~masked_variance_pseudo_psu, weights = ~WTMEC4YR_norm,
            strata = ~masked_variance_pseudo_stratum, data=., nest=TRUE)

# survey weighted table
df_svy %>%
  tbl_svysummary(
    by = gender,
    include = -c(SEQN, data_release_cycle, masked_variance_pseudo_psu, masked_variance_pseudo_stratum, full_sample_2_year_interview_weight, full_sample_2_year_mec_exam_weight, WTMEC4YR, WTMEC4YR_norm),
    statistic = list(
      all_continuous() ~ "{mean} ({sd})",
      all_categorical() ~ "{n} ({p}%)"
    ),
    digits = all_continuous() ~ 1,
    missing_text = "Missing",
  ) %>%
  add_overall() %>%
  add_p()

Characteristic	Overall N = 8,968 1	Female N = 4,699 1	Male N = 4,269 1	p-value 2
Wear minutes	929.1 (123.6)	916.2 (119.5)	943.3 (126.4)	<0.001
Age	45.9 (17.3)	46.4 (17.6)	45.3 (16.9)	0.022
Total AC	260,231.6 (183,008.2)	231,075.7 (200,451.9)	292,324.9 (155,469.4)	<0.001
Mean AC	180.7 (127.1)	160.5 (139.2)	203.0 (108.0)	<0.001
MVPA minutes	22.7 (28.3)	16.9 (29.3)	29.1 (25.8)	<0.001
Moderate minutes	21.7 (22.4)	16.0 (18.8)	27.9 (24.4)	<0.001
Vigorous minutes	1.0 (10.6)	0.9 (14.1)	1.2 (4.2)	<0.001
MVPA (10 min bout)	4.2 (19.5)	3.8 (25.3)	4.6 (9.7)	<0.001
Vig (10 min bout)	0.5 (9.6)	0.5 (12.9)	0.4 (2.9)	<0.001
Time <5	895.3 (128.2)	899.8 (124.0)	890.4 (132.6)	0.028
Time [5,10]	28.7 (9.6)	29.5 (9.6)	27.7 (9.5)	<0.001
Time (10,25]	51.8 (16.1)	53.2 (16.2)	50.3 (16.0)	<0.001
Time (25,50]	52.3 (15.5)	53.3 (15.4)	51.1 (15.5)	<0.001
Time (50,75]	36.7 (10.9)	37.3 (10.7)	35.9 (11.0)	<0.001
Time (75,100]	28.8 (8.6)	29.4 (8.4)	28.1 (8.7)	<0.001
Time (100,150]	44.1 (13.7)	45.1 (13.0)	42.9 (14.4)	<0.001
Time (150,200]	33.2 (10.3)	34.2 (10.2)	32.1 (10.3)	<0.001
Time (200,500]	113.2 (36.9)	116.5 (37.4)	109.5 (36.0)	<0.001
Time (500,1000]	80.5 (35.5)	79.1 (35.1)	82.0 (35.9)	0.012
Time (1000,2000]	52.2 (33.5)	45.0 (28.3)	60.2 (36.9)	<0.001
Time (2000,4000]	19.1 (18.1)	14.2 (13.4)	24.5 (20.9)	<0.001
Time (4000,6000]	3.1 (9.7)	2.2 (11.7)	4.0 (6.7)	<0.001
Time >6000	1.0 (9.9)	0.8 (13.0)	1.2 (4.2)	<0.001
Missing AC	502.3 (133.2)	501.0 (141.3)	517.6 (3.4)
Missing	8,959	4,691	4,268
Median AC	5.8 (88.8)	5.7 (121.7)	6.0 (16.7)	0.001
log10 mean AC	0.9 (0.2)	0.8 (0.2)	0.9 (0.2)	<0.001
log10 median AC	0.2 (0.4)	0.2 (0.3)	0.3 (0.4)	0.007
log10 total AC	1,247.6 (336.2)	1,220.2 (319.1)	1,277.8 (351.6)	<0.001
1 Mean (SD)
2 Design-based KruskalWallis test

Medians

Code

# unweighted, use median
df_small %>% 
  tbl_summary(
    by = gender,
    include = -c(SEQN, data_release_cycle, masked_variance_pseudo_psu, masked_variance_pseudo_stratum, full_sample_2_year_interview_weight, full_sample_2_year_mec_exam_weight),
    statistic = list(
      all_continuous() ~ "{median} ({IQR})",
      all_categorical() ~ "{n} ({p}%)"
    ),
    digits = all_continuous() ~ 1,
    missing_text = "Missing",
  ) %>%
  add_overall() %>%
  add_p()

Characteristic	Overall N = 8,968 1	Female N = 4,671 1	Male N = 4,297 1	p-value 2
Wear minutes	920.4 (158.3)	904.0 (155.7)	939.6 (157.1)	<0.001
Age	45.0 (35.0)	44.0 (34.0)	46.0 (35.0)	<0.001
Total AC	225,022.3 (166,694.6)	203,857.9 (131,907.2)	259,047.7 (204,421.5)	<0.001
Mean AC	156.3 (115.8)	141.6 (91.6)	179.9 (142.0)	<0.001
MVPA minutes	13.9 (24.9)	9.9 (17.3)	20.8 (33.5)	<0.001
Moderate minutes	13.6 (24.0)	9.7 (16.8)	20.3 (31.9)	<0.001
Vigorous minutes	0.0 (0.3)	0.0 (0.1)	0.0 (0.5)	<0.001
MVPA (10 min bout)	0.0 (3.3)	0.0 (1.8)	0.2 (5.0)	<0.001
Vig (10 min bout)	0.0 (0.0)	0.0 (0.0)	0.0 (0.0)	<0.001
Time <5	902.8 (179.0)	909.0 (171.4)	894.9 (190.3)	<0.001
Time [5,10]	27.4 (12.7)	28.3 (12.9)	26.5 (12.3)	<0.001
Time (10,25]	49.7 (20.9)	50.7 (21.2)	48.6 (20.7)	<0.001
Time (25,50]	50.4 (19.9)	51.3 (20.1)	49.3 (20.0)	<0.001
Time (50,75]	35.4 (14.1)	36.0 (14.1)	34.7 (14.2)	<0.001
Time (75,100]	28.0 (11.0)	28.6 (11.0)	27.3 (11.0)	<0.001
Time (100,150]	43.0 (17.2)	44.0 (17.5)	42.0 (16.9)	<0.001
Time (150,200]	32.7 (13.6)	33.6 (13.6)	31.7 (13.4)	<0.001
Time (200,500]	110.9 (48.4)	113.7 (49.8)	107.4 (46.3)	<0.001
Time (500,1000]	74.6 (48.2)	73.3 (46.5)	76.5 (50.7)	<0.001
Time (1000,2000]	42.7 (40.9)	37.5 (34.4)	49.8 (49.8)	<0.001
Time (2000,4000]	12.7 (20.8)	9.4 (15.3)	18.2 (26.4)	<0.001
Time (4000,6000]	0.5 (2.7)	0.3 (1.3)	1.0 (4.5)	<0.001
Time >6000	0.0 (0.3)	0.0 (0.1)	0.0 (0.5)	<0.001
Missing AC	529.0 (41.5)	544.0 (45.8)	519.0 (3.5)	0.3
Missing	8,957	4,663	4,294
Median AC	0.1 (3.4)	0.0 (2.4)	0.3 (4.9)	<0.001
log10 mean AC	0.8 (0.3)	0.8 (0.3)	0.9 (0.4)	<0.001
log10 median AC	0.0 (0.3)	0.0 (0.3)	0.1 (0.4)	<0.001
log10 total AC	1,216.0 (460.9)	1,188.1 (420.0)	1,253.0 (506.0)	<0.001
1 Median (IQR)
2 Wilcoxon rank sum test

Stratified by wave

2003-2004

Code

df_svy =
  df_small %>% 
  filter(data_release_cycle == 3) %>%
  mutate(WTMEC4YR = full_sample_2_year_mec_exam_weight,
         WTMEC4YR_norm = WTMEC4YR/mean(WTMEC4YR, na.rm = TRUE)) %>%
  svydesign(ids = ~masked_variance_pseudo_psu, weights = ~WTMEC4YR_norm,
            strata = ~masked_variance_pseudo_stratum, data=., nest=TRUE)

# survey weighted table
df_svy %>%
  tbl_svysummary(
    by = gender,
    include = -c(SEQN, data_release_cycle, masked_variance_pseudo_psu, masked_variance_pseudo_stratum, full_sample_2_year_interview_weight, full_sample_2_year_mec_exam_weight, WTMEC4YR, WTMEC4YR_norm),
    statistic = list(
      all_continuous() ~ "{mean} ({sd})",
      all_categorical() ~ "{n} ({p}%)"
    ),
    digits = all_continuous() ~ 0,
    missing_text = "Missing",
  ) %>%
  add_overall() %>%
  add_p()

Characteristic	Overall N = 4,503 1	Female N = 2,349 1	Male N = 2,154 1	p-value 2
Wear minutes	926 (121)	915 (118)	938 (124)	<0.001
Age	46 (17)	46 (18)	45 (17)	0.069
Total AC	261,985 (222,063)	233,780 (262,098)	292,741 (162,426)	<0.001
Mean AC	182 (154)	162 (182)	203 (113)	<0.001
MVPA minutes	23 (33)	17 (38)	29 (26)	<0.001
Moderate minutes	22 (24)	16 (22)	28 (24)	<0.001
Vigorous minutes	1 (15)	1 (20)	1 (5)	<0.001
MVPA (10 min bout)	5 (26)	4 (35)	5 (10)	<0.001
Vig (10 min bout)	1 (13)	1 (18)	0 (3)	<0.001
Time <5	896 (127)	898 (124)	894 (131)	0.4
Time [5,10]	28 (10)	30 (10)	27 (9)	<0.001
Time (10,25]	51 (16)	53 (16)	49 (16)	<0.001
Time (25,50]	52 (16)	53 (16)	50 (15)	<0.001
Time (50,75]	36 (11)	37 (11)	35 (11)	<0.001
Time (75,100]	29 (9)	30 (9)	28 (9)	<0.001
Time (100,150]	44 (14)	46 (13)	43 (15)	<0.001
Time (150,200]	33 (10)	35 (10)	32 (10)	<0.001
Time (200,500]	114 (37)	118 (37)	110 (35)	<0.001
Time (500,1000]	80 (35)	79 (35)	82 (35)	0.038
Time (1000,2000]	52 (34)	44 (28)	60 (37)	<0.001
Time (2000,4000]	19 (18)	14 (14)	25 (21)	<0.001
Time (4000,6000]	3 (13)	2 (16)	4 (7)	<0.001
Time >6000	1 (14)	1 (18)	1 (5)	<0.001
Missing AC
95	1 (24%)	1 (28%)	0 (0%)
512	0 (4.6%)	0 (0%)	0 (29%)
519	0 (11%)	0 (0%)	0 (71%)
529	1 (43%)	1 (51%)	0 (0%)
563	1 (18%)	1 (21%)	0 (0%)
Missing	4,500	2,346	2,154
Median AC	7 (125)	8 (172)	6 (13)	0.11
log10 mean AC	1 (0)	1 (0)	1 (0)	0.007
log10 median AC	0 (0)	0 (0)	0 (0)	0.2
log10 total AC	1,248 (335)	1,225 (321)	1,272 (348)	0.007
1 Mean (SD); n (%)
2 Design-based KruskalWallis test

2005-2006

Code

df_svy =
  df_small %>% 
  filter(data_release_cycle == 4) %>%
  mutate(WTMEC4YR = full_sample_2_year_mec_exam_weight,
         WTMEC4YR_norm = WTMEC4YR/mean(WTMEC4YR, na.rm = TRUE)) %>%
  svydesign(ids = ~masked_variance_pseudo_psu, weights = ~WTMEC4YR_norm,
            strata = ~masked_variance_pseudo_stratum, data=., nest=TRUE)

# survey weighted table
df_svy %>%
  tbl_svysummary(
    by = gender,
    include = -c(SEQN, data_release_cycle, masked_variance_pseudo_psu, masked_variance_pseudo_stratum, full_sample_2_year_interview_weight, full_sample_2_year_mec_exam_weight, WTMEC4YR, WTMEC4YR_norm),
    statistic = list(
      all_continuous() ~ "{mean} ({sd})",
      all_categorical() ~ "{n} ({p}%)"
    ),
    digits = all_continuous() ~ 0,
    missing_text = "Missing",
  ) %>%
  add_overall() %>%
  add_p()

Characteristic	Overall N = 4,465 1	Female N = 2,350 1	Male N = 2,115 1	p-value 2
Wear minutes	932 (126)	917 (121)	949 (129)	<0.001
Age	46 (17)	47 (18)	46 (17)	0.2
Total AC	258,490 (133,323)	228,413 (110,014)	291,907 (148,198)	<0.001
Mean AC	180 (93)	159 (76)	203 (103)	<0.001
MVPA minutes	22 (22)	16 (16)	29 (26)	<0.001
Moderate minutes	21 (21)	16 (15)	28 (24)	<0.001
Vigorous minutes	1 (3)	1 (2)	1 (4)	<0.001
MVPA (10 min bout)	4 (8)	3 (7)	4 (9)	<0.001
Vig (10 min bout)	0 (2)	0 (2)	0 (2)	0.038
Time <5	895 (129)	902 (124)	887 (134)	0.021
Time [5,10]	29 (10)	30 (9)	28 (10)	<0.001
Time (10,25]	52 (16)	53 (16)	51 (16)	0.001
Time (25,50]	53 (15)	53 (15)	52 (16)	0.043
Time (50,75]	37 (11)	37 (11)	37 (11)	0.046
Time (75,100]	29 (9)	29 (8)	28 (9)	0.005
Time (100,150]	44 (13)	45 (13)	43 (13)	0.003
Time (150,200]	33 (10)	34 (10)	32 (10)	<0.001
Time (200,500]	112 (37)	115 (37)	109 (37)	<0.001
Time (500,1000]	81 (36)	79 (36)	82 (36)	0.13
Time (1000,2000]	53 (34)	46 (29)	61 (36)	<0.001
Time (2000,4000]	19 (18)	14 (13)	24 (21)	<0.001
Time (4000,6000]	3 (6)	2 (5)	4 (6)	<0.001
Time >6000	1 (3)	1 (2)	1 (4)	<0.001
Missing AC
501	3 (46%)	3 (47%)	0 (0%)
519	0 (3.6%)	0 (0%)	0 (100%)
537	1 (12%)	1 (12%)	0 (0%)
551	1 (18%)	1 (19%)	0 (0%)
582	1 (11%)	1 (12%)	0 (0%)
662	1 (9.5%)	1 (9.8%)	0 (0%)
Missing	4,459	2,344	2,115
Median AC	5 (15)	4 (10)	6 (20)	0.003
log10 mean AC	1 (0)	1 (0)	1 (0)	<0.001
log10 median AC	0 (0)	0 (0)	0 (0)	0.006
log10 total AC	1,248 (337)	1,215 (317)	1,284 (355)	<0.001
1 Mean (SD); n (%)
2 Design-based KruskalWallis test

Survey-Weighted Physical Activity Summaries by Sex and Age

Code

survey_design =
  df_small %>%
  mutate(age_cat = cut(age_in_years_at_screening,
                       breaks=c(18, 30, seq(40, 80, 10)))) %>%
  mutate(weight = full_sample_2_year_mec_exam_weight / 2,
         weight_norm = weight / mean(weight)) %>%
  ungroup()  %>%
  survey::svydesign(
    id = ~ masked_variance_pseudo_psu,
    strata = ~ masked_variance_pseudo_stratum,
    weights = ~ weight_norm,
    data = .,
    nest = TRUE
  )

calc_by_age =
  function(var, df) {
    formula = as.formula(paste("~", var))
    total_by_age_gender = svyby(formula,
                                ~ age_cat + gender,
                                df,
                                svymean) %>%
      rename(mean = contains(var)) %>%
      mutate(metric = var)
  }

means_df =
  map_dfr(.x = pa_old %>% select(contains("bout"), starts_with("x"), total, mean, mvpa, mod, vig) %>% colnames(),
          .f = calc_by_age, df = survey_design)


small = 
  means_df %>%
  pivot_wider(names_from = gender, values_from = c(mean, se)) %>% 
  mutate(val_female = paste0(round(mean_Female), " (", round(mean_Female - 1.96 * se_Female), ",", round(mean_Female + 1.96 * se_Female), ")"),
         val_male = paste0(round(mean_Male), " (", round(mean_Male - 1.96 * se_Female), ",", round(mean_Male + 1.96 * se_Male), ")"))

# test = svyttest(mvpa_montoye ~ gender, survey_design)
pvals = c()
for(i in 1:nrow(small)){
  row = means_df[i,]
  cat = row$age_cat
  var = row$metric
  survey_design =
    df_small %>%
    mutate(age_cat = cut(age_in_years_at_screening,
                         breaks=c(18, 30, seq(40, 80, 10)))) %>%
    filter(age_cat == cat) %>%
    mutate(weight = full_sample_2_year_mec_exam_weight / 2,
           weight_norm = weight / mean(weight)) %>%
    ungroup()  %>%
    survey::svydesign(
      id = ~ masked_variance_pseudo_psu,
      strata = ~ masked_variance_pseudo_stratum,
      weights = ~ weight_norm,
      data = .,
      nest = TRUE
    )
  formula = as.formula(paste(var,  "~ gender"))
  test = svyttest(formula, survey_design)
  pvals = append(pvals, unname(test$p.value))
}

small %>% 
  select(metric, age_cat, Female = val_female, Male = val_male, mean_Female, mean_Male) %>% 
  bind_cols(p = pvals) %>% 
  mutate(pval = if_else(p < 0.001, "<0.001", format.pval(p, digits = 2))) %>% 
  # mutate(metric = factor(metric, 
  #                        labels = c("MVPA (min)", 
  #                                   "VPA (min)", "Time <1000",
  #                                   "Time [1000,2000]", "Time (2000,3000]", "Time (3000,4000]", "Time (4000,5000]", "Time (5000,6000]", "Time (6000,7000]", "Time (7000,8000]", "Time (8000,9000]", "Time >9000", "Total AC", "Mean AC", "Total MIMS"))) %>% 
  group_by(metric) %>% 
  gt() %>% 
  cols_hide(c("mean_Female", "mean_Male", "p")) %>% 
  tab_style(
    # Style to apply when conditions are met
    style = cell_fill(color = "lightpink"),
    locations = cells_body(
      rows = mean_Female > mean_Male & p < 0.05
    )
  ) %>% 
  tab_style(
    # Style to apply when conditions are met
    style = cell_fill(color = "lightblue"),
    locations = cells_body(
      rows = mean_Female < mean_Male & p < 0.05
    )
  )

age_cat	Female	Male	pval
mvpa_bout
(18,30]	4 (3,4)	6 (5,7)	<0.001
(30,40]	4 (3,4)	5 (4,5)	0.01738
(40,50]	4 (3,5)	5 (4,6)	0.27386
(50,60]	3 (2,4)	4 (3,5)	0.03474
(60,70]	3 (2,4)	4 (3,5)	0.38607
(70,80]	7 (-3,18)	2 (-8,3)	0.32701
vig_bout
(18,30]	0 (0,1)	0 (0,1)	<0.001
(30,40]	0 (0,1)	1 (0,1)	0.01738
(40,50]	0 (0,0)	1 (1,1)	0.27386
(50,60]	0 (0,0)	0 (0,1)	0.03474
(60,70]	0 (0,0)	0 (0,0)	0.38607
(70,80]	3 (-3,8)	0 (-5,0)	0.32701
x0_5
(18,30]	914 (904,923)	877 (868,888)	0.69216
(30,40]	879 (871,887)	855 (847,868)	0.49795
(40,50]	871 (858,885)	870 (857,885)	0.03857
(50,60]	878 (868,887)	886 (876,903)	0.04307
(60,70]	916 (903,928)	934 (922,949)	0.21096
(70,80]	949 (931,966)	984 (967,997)	0.32759
x5_10
(18,30]	28 (27,28)	26 (25,27)	0.69216
(30,40]	28 (27,29)	26 (25,27)	0.49795
(40,50]	30 (29,31)	28 (27,29)	0.03857
(50,60]	31 (30,32)	29 (28,30)	0.04307
(60,70]	31 (30,32)	30 (29,30)	0.21096
(70,80]	32 (31,33)	30 (28,31)	0.32759
x10_25
(18,30]	50 (49,52)	48 (47,49)	<0.001
(30,40]	52 (50,54)	49 (48,51)	0.00583
(40,50]	54 (52,56)	51 (49,53)	0.93900
(50,60]	55 (53,57)	51 (49,54)	0.31344
(60,70]	53 (52,55)	52 (50,54)	0.00617
(70,80]	56 (53,58)	52 (49,53)	<0.001
x25_50
(18,30]	51 (50,52)	50 (49,51)	<0.001
(30,40]	53 (51,55)	51 (50,53)	0.00583
(40,50]	54 (53,56)	52 (51,54)	0.93900
(50,60]	55 (53,56)	52 (50,54)	0.31344
(60,70]	53 (52,55)	51 (49,53)	0.00617
(70,80]	54 (52,56)	51 (49,52)	<0.001
x50_75
(18,30]	36 (35,36)	35 (35,36)	0.00273
(30,40]	38 (37,39)	36 (35,38)	0.01111
(40,50]	38 (37,39)	37 (36,38)	0.00889
(50,60]	39 (38,40)	37 (36,38)	0.02781
(60,70]	37 (36,38)	35 (34,36)	0.04774
(70,80]	37 (36,39)	35 (34,36)	0.02226
x75_100
(18,30]	28 (27,28)	28 (27,29)	0.00273
(30,40]	30 (29,31)	28 (27,29)	0.01111
(40,50]	30 (29,31)	29 (28,30)	0.00889
(50,60]	31 (30,31)	29 (28,30)	0.02781
(60,70]	29 (28,30)	27 (26,28)	0.04774
(70,80]	29 (28,30)	27 (26,27)	0.02226
x100_150
(18,30]	43 (42,44)	43 (42,44)	0.00216
(30,40]	46 (44,47)	43 (42,45)	0.02779
(40,50]	46 (44,48)	44 (43,46)	0.01928
(50,60]	47 (46,48)	44 (43,45)	0.01127
(60,70]	45 (44,47)	41 (40,43)	0.07634
(70,80]	45 (44,47)	40 (38,41)	0.00789
x150_200
(18,30]	32 (31,33)	32 (31,33)	0.00216
(30,40]	35 (34,35)	33 (32,34)	0.02779
(40,50]	35 (34,36)	32 (31,34)	0.01928
(50,60]	36 (35,37)	33 (32,34)	0.01127
(60,70]	35 (34,36)	31 (30,32)	0.07634
(70,80]	34 (33,36)	30 (28,31)	0.00789
x200_500
(18,30]	109 (106,112)	110 (107,114)	0.07747
(30,40]	119 (116,121)	114 (111,118)	0.12478
(40,50]	121 (118,125)	111 (107,114)	0.10130
(50,60]	125 (121,128)	113 (110,117)	0.01986
(60,70]	121 (118,124)	107 (103,110)	0.01089
(70,80]	111 (106,117)	97 (91,101)	0.00596
x500_1e_03
(18,30]	79 (76,82)	85 (82,89)	0.07747
(30,40]	87 (83,90)	92 (88,96)	0.12478
(40,50]	88 (84,91)	86 (83,89)	0.10130
(50,60]	85 (81,88)	84 (81,88)	0.01986
(60,70]	74 (71,77)	73 (70,77)	0.01089
(70,80]	57 (53,61)	58 (54,62)	0.00596
x1e_03_2e_03
(18,30]	50 (48,52)	67 (65,71)	0.54919
(30,40]	54 (51,57)	74 (71,78)	0.16870
(40,50]	52 (50,55)	66 (64,70)	0.14382
(50,60]	45 (43,48)	59 (56,62)	0.00387
(60,70]	35 (32,37)	43 (41,47)	0.00318
(70,80]	22 (20,24)	28 (26,31)	0.00730
x2e_03_4e_03
(18,30]	18 (16,19)	31 (30,34)	0.54919
(30,40]	17 (16,18)	31 (30,34)	0.16870
(40,50]	17 (16,18)	27 (26,30)	0.14382
(50,60]	13 (11,14)	20 (19,22)	0.00387
(60,70]	9 (8,10)	14 (13,15)	0.00318
(70,80]	6 (5,7)	8 (7,9)	0.00730
x4e_03_6e_03
(18,30]	3 (2,3)	6 (5,7)	0.58355
(30,40]	3 (2,3)	5 (4,5)	0.06021
(40,50]	2 (2,3)	4 (3,5)	0.15626
(50,60]	2 (1,2)	3 (3,4)	0.00150
(60,70]	1 (1,1)	2 (2,3)	<0.001
(70,80]	3 (-2,8)	1 (-4,1)	<0.001
x6e_03_inf
(18,30]	1 (1,1)	2 (1,2)	0.58355
(30,40]	1 (1,1)	1 (1,2)	0.06021
(40,50]	1 (1,1)	1 (1,2)	0.15626
(50,60]	0 (0,0)	1 (1,1)	0.00150
(60,70]	0 (0,0)	0 (0,0)	<0.001
(70,80]	3 (-3,8)	0 (-5,0)	<0.001
total
(18,30]	249414 (240645,258183)	338721 (329952,352897)	0.87103
(30,40]	259369 (251096,267642)	344614 (336341,357723)	0.03300
(40,50]	256144 (247892,264396)	313436 (305184,328690)	0.13275
(50,60]	227431 (216765,238097)	273748 (263082,288199)	0.00158
(60,70]	188786 (179801,197772)	213610 (204624,223538)	<0.001
(70,80]	182335 (109907,254763)	154885 (82457,163860)	<0.001
mean
(18,30]	173 (167,179)	235 (229,245)	0.87103
(30,40]	180 (174,186)	239 (234,248)	0.03300
(40,50]	178 (172,184)	218 (212,228)	0.13275
(50,60]	158 (151,165)	190 (183,200)	0.00158
(60,70]	131 (125,137)	148 (142,155)	<0.001
(70,80]	127 (76,177)	108 (57,114)	<0.001
mvpa
(18,30]	21 (19,23)	38 (37,42)	0.72753
(30,40]	20 (18,21)	37 (35,39)	0.01495
(40,50]	19 (18,21)	32 (31,35)	0.01013
(50,60]	14 (13,16)	24 (22,26)	<0.001
(60,70]	10 (9,12)	16 (14,17)	<0.001
(70,80]	12 (1,22)	9 (-2,10)	<0.001
mod
(18,30]	20 (18,22)	37 (35,40)	0.72753
(30,40]	19 (18,20)	35 (34,38)	0.01495
(40,50]	19 (17,20)	31 (29,33)	0.01013
(50,60]	14 (12,16)	23 (21,25)	<0.001
(60,70]	10 (9,11)	15 (14,17)	<0.001
(70,80]	9 (4,13)	8 (4,10)	<0.001
vig
(18,30]	1 (1,1)	2 (1,2)	0.67430
(30,40]	1 (1,1)	1 (1,2)	0.05279
(40,50]	1 (1,1)	1 (1,2)	0.00163
(50,60]	0 (0,0)	1 (1,1)	<0.001
(60,70]	0 (0,0)	0 (0,0)	<0.001
(70,80]	3 (-3,9)	0 (-6,0)	<0.001

Medians

Code

survey_design =
  df_small %>%
  mutate(age_cat = cut(age_in_years_at_screening,
                       breaks=c(18, 30, seq(40, 80, 10)))) %>% 
  filter(!is.na(age_cat))
calc_by_age =
  function(var, df) {
    df %>% 
      select(v = all_of(var), age_cat, gender) %>% 
      group_by(age_cat, gender) %>% 
      summarize(across(v, ~median(.x, na.rm = TRUE))) %>% 
      rename(median = contains("v")) %>%
      mutate(metric = var)
  }

means_df =
  map_dfr(.x = pa_old %>% select(contains("bout"), starts_with("x"), total, mean, mvpa, mod, vig) %>% colnames(),
          .f = calc_by_age, df = survey_design)


small = 
  means_df %>%
  pivot_wider(names_from = gender, values_from = c(median)) 

small %>% 
  select(metric, age_cat, Female, Male) %>% 
  group_by(metric) %>% 
  gt() %>% 
  tab_style(
    # Style to apply when conditions are met
    style = cell_fill(color = "lightpink"),
    locations = cells_body(
      rows = Female > Male 
    )
  ) %>% 
  tab_style(
    # Style to apply when conditions are met
    style = cell_fill(color = "lightblue"),
    locations = cells_body(
      rows = Female < Male 
    )
  )

age_cat	Female	Male
mvpa_bout
(18,30]	0.000000e+00	2.000000e+00
(30,40]	0.000000e+00	1.000000e+00
(40,50]	0.000000e+00	5.357143e-01
(50,60]	0.000000e+00	0.000000e+00
(60,70]	0.000000e+00	0.000000e+00
(70,80]	0.000000e+00	0.000000e+00
vig_bout
(18,30]	0.000000e+00	0.000000e+00
(30,40]	0.000000e+00	0.000000e+00
(40,50]	0.000000e+00	0.000000e+00
(50,60]	0.000000e+00	0.000000e+00
(60,70]	0.000000e+00	0.000000e+00
(70,80]	0.000000e+00	0.000000e+00
x0_5
(18,30]	9.228571e+02	8.660000e+02
(30,40]	8.790000e+02	8.396000e+02
(40,50]	8.710000e+02	8.587143e+02
(50,60]	8.750000e+02	8.814167e+02
(60,70]	9.082857e+02	9.322857e+02
(70,80]	9.570000e+02	9.750000e+02
x5_10
(18,30]	2.620000e+01	2.500000e+01
(30,40]	2.680000e+01	2.500000e+01
(40,50]	2.860000e+01	2.600000e+01
(50,60]	2.900000e+01	2.750000e+01
(60,70]	2.973214e+01	2.781667e+01
(70,80]	3.166667e+01	2.824286e+01
x10_25
(18,30]	4.783333e+01	4.633333e+01
(30,40]	5.016667e+01	4.757143e+01
(40,50]	5.142857e+01	4.942857e+01
(50,60]	5.314286e+01	4.966667e+01
(60,70]	5.171429e+01	4.924286e+01
(70,80]	5.442857e+01	4.992857e+01
x25_50
(18,30]	4.820000e+01	4.814286e+01
(30,40]	5.142857e+01	5.000000e+01
(40,50]	5.320000e+01	5.080000e+01
(50,60]	5.400000e+01	5.100000e+01
(60,70]	5.214286e+01	4.824286e+01
(70,80]	5.342857e+01	4.942857e+01
x50_75
(18,30]	3.420000e+01	3.420000e+01
(30,40]	3.633333e+01	3.600000e+01
(40,50]	3.680000e+01	3.558571e+01
(50,60]	3.800000e+01	3.571429e+01
(60,70]	3.714286e+01	3.330952e+01
(70,80]	3.742857e+01	3.400000e+01
x75_100
(18,30]	2.733333e+01	2.760000e+01
(30,40]	2.885714e+01	2.814286e+01
(40,50]	2.966667e+01	2.800000e+01
(50,60]	2.983333e+01	2.816667e+01
(60,70]	2.930952e+01	2.561905e+01
(70,80]	2.933333e+01	2.630952e+01
x100_150
(18,30]	4.171429e+01	4.200000e+01
(30,40]	4.471429e+01	4.328571e+01
(40,50]	4.542857e+01	4.333333e+01
(50,60]	4.614286e+01	4.315476e+01
(60,70]	4.638095e+01	4.000000e+01
(70,80]	4.557143e+01	4.042857e+01
x150_200
(18,30]	3.120000e+01	3.200000e+01
(30,40]	3.440000e+01	3.300000e+01
(40,50]	3.471429e+01	3.253571e+01
(50,60]	3.614286e+01	3.284524e+01
(60,70]	3.550000e+01	3.057143e+01
(70,80]	3.500000e+01	3.000000e+01
x200_500
(18,30]	1.090000e+02	1.094286e+02
(30,40]	1.160000e+02	1.155000e+02
(40,50]	1.217500e+02	1.121667e+02
(50,60]	1.240000e+02	1.130000e+02
(60,70]	1.200714e+02	1.044286e+02
(70,80]	1.104286e+02	9.821429e+01
x500_1e_03
(18,30]	7.314286e+01	8.233333e+01
(30,40]	8.300000e+01	9.100000e+01
(40,50]	8.483333e+01	8.561905e+01
(50,60]	8.271429e+01	8.161905e+01
(60,70]	6.960000e+01	7.138095e+01
(70,80]	4.866667e+01	5.714286e+01
x1e_03_2e_03
(18,30]	4.250000e+01	5.880000e+01
(30,40]	4.800000e+01	6.750000e+01
(40,50]	4.614286e+01	6.021429e+01
(50,60]	4.050000e+01	5.185714e+01
(60,70]	2.771429e+01	3.775714e+01
(70,80]	1.471429e+01	2.230952e+01
x2e_03_4e_03
(18,30]	1.342857e+01	2.833333e+01
(30,40]	1.300000e+01	2.633333e+01
(40,50]	1.271429e+01	2.300000e+01
(50,60]	8.285714e+00	1.742857e+01
(60,70]	4.714286e+00	9.000000e+00
(70,80]	2.166667e+00	4.000000e+00
x4e_03_6e_03
(18,30]	7.142857e-01	3.500000e+00
(30,40]	5.714286e-01	2.666667e+00
(40,50]	4.285714e-01	1.535714e+00
(50,60]	2.857143e-01	6.904762e-01
(60,70]	1.428571e-01	2.500000e-01
(70,80]	0.000000e+00	7.142857e-02
x6e_03_inf
(18,30]	0.000000e+00	4.285714e-01
(30,40]	0.000000e+00	1.666667e-01
(40,50]	0.000000e+00	0.000000e+00
(50,60]	0.000000e+00	0.000000e+00
(60,70]	0.000000e+00	0.000000e+00
(70,80]	0.000000e+00	0.000000e+00
total
(18,30]	2.209614e+05	3.241940e+05
(30,40]	2.389592e+05	3.258449e+05
(40,50]	2.353852e+05	3.032826e+05
(50,60]	2.145245e+05	2.588716e+05
(60,70]	1.707527e+05	2.015814e+05
(70,80]	1.222527e+05	1.412056e+05
mean
(18,30]	1.534454e+02	2.251347e+02
(30,40]	1.659439e+02	2.262812e+02
(40,50]	1.634619e+02	2.106129e+02
(50,60]	1.489753e+02	1.797719e+02
(60,70]	1.185783e+02	1.399871e+02
(70,80]	8.489772e+01	9.805944e+01
mvpa
(18,30]	1.442857e+01	3.500000e+01
(30,40]	1.366667e+01	3.116667e+01
(40,50]	1.400000e+01	2.657143e+01
(50,60]	8.400000e+00	1.942857e+01
(60,70]	4.732143e+00	9.571429e+00
(70,80]	2.285714e+00	3.857143e+00
mod
(18,30]	1.400000e+01	3.316667e+01
(30,40]	1.357143e+01	3.000000e+01
(40,50]	1.360000e+01	2.564286e+01
(50,60]	8.400000e+00	1.921429e+01
(60,70]	4.714286e+00	9.500000e+00
(70,80]	2.285714e+00	3.803571e+00
vig
(18,30]	0.000000e+00	4.285714e-01
(30,40]	0.000000e+00	1.666667e-01
(40,50]	0.000000e+00	0.000000e+00
(50,60]	0.000000e+00	0.000000e+00
(60,70]	0.000000e+00	0.000000e+00
(70,80]	0.000000e+00	0.000000e+00

Stratified by wave

2003-2004

Code

survey_design =
  df_small %>%
  filter(data_release_cycle == 3) %>%
  mutate(age_cat = cut(age_in_years_at_screening,
                       breaks=c(18, 30, seq(40, 80, 10)))) %>%
  mutate(weight = full_sample_2_year_mec_exam_weight,
         weight_norm = weight / mean(weight)) %>%
  ungroup()  %>%
  survey::svydesign(
    id = ~ masked_variance_pseudo_psu,
    strata = ~ masked_variance_pseudo_stratum,
    weights = ~ weight_norm,
    data = .,
    nest = TRUE
  )

calc_by_age =
  function(var, df) {
    formula = as.formula(paste("~", var))
    total_by_age_gender = svyby(formula,
                                ~ age_cat + gender,
                                df,
                                svymean) %>%
      rename(mean = contains(var)) %>%
      mutate(metric = var)
  }

means_df =
  map_dfr(.x = pa_old %>% select(contains("bout"), starts_with("x"), total, mean, mvpa, mod, vig) %>% colnames(),
          .f = calc_by_age, df = survey_design)


small = 
  means_df %>%
  pivot_wider(names_from = gender, values_from = c(mean, se)) %>% 
  mutate(val_female = paste0(round(mean_Female), " (", round(mean_Female - 1.96 * se_Female), ",", round(mean_Female + 1.96 * se_Female), ")"),
         val_male = paste0(round(mean_Male), " (", round(mean_Male - 1.96 * se_Female), ",", round(mean_Male + 1.96 * se_Male), ")"))

# test = svyttest(mvpa_montoye ~ gender, survey_design)
pvals = c()
for(i in 1:nrow(small)){
  row = means_df[i,]
  cat = row$age_cat
  var = row$metric
  survey_design =
    df_small %>%
    filter(age_in_years_at_screening >= 18 & data_release_cycle == 3) %>%
    mutate(age_cat = cut(age_in_years_at_screening,
                         breaks=c(18, 30, seq(40, 80, 10)))) %>%
    filter(age_cat == cat) %>%
    mutate(weight = full_sample_2_year_mec_exam_weight / 2,
           weight_norm = weight / mean(weight)) %>%
    ungroup()  %>%
    survey::svydesign(
      id = ~ masked_variance_pseudo_psu,
      strata = ~ masked_variance_pseudo_stratum,
      weights = ~ weight_norm,
      data = .,
      nest = TRUE
    )
  formula = as.formula(paste(var,  "~ gender"))
  test = svyttest(formula, survey_design)
  pvals = append(pvals, unname(test$p.value))
}

small %>% 
  select(metric, age_cat, Female = val_female, Male = val_male, mean_Female, mean_Male) %>% 
  bind_cols(p = pvals) %>% 
  mutate(pval = if_else(p < 0.001, "<0.001", format.pval(p, digits = 2))) %>% 
  # mutate(metric = factor(metric, 
  #                        labels = c("MVPA (min)", 
  #                                   "VPA (min)", "Time <1000",
  #                                   "Time [1000,2000]", "Time (2000,3000]", "Time (3000,4000]", "Time (4000,5000]", "Time (5000,6000]", "Time (6000,7000]", "Time (7000,8000]", "Time (8000,9000]", "Time >9000", "Total AC", "Mean AC", "Total MIMS"))) %>% 
  group_by(metric) %>% 
  gt() %>% 
  cols_hide(c("mean_Female", "mean_Male", "p")) %>% 
  tab_style(
    # Style to apply when conditions are met
    style = cell_fill(color = "lightpink"),
    locations = cells_body(
      rows = mean_Female > mean_Male & p < 0.05
    )
  ) %>% 
  tab_style(
    # Style to apply when conditions are met
    style = cell_fill(color = "lightblue"),
    locations = cells_body(
      rows = mean_Female < mean_Male & p < 0.05
    )
  )

age_cat	Female	Male	pval
mvpa_bout
(18,30]	4 (3,6)	6 (4,7)	0.02386
(30,40]	3 (2,4)	5 (4,6)	<0.001
(40,50]	4 (3,5)	5 (4,7)	0.19637
(50,60]	3 (2,4)	5 (3,6)	0.07122
(60,70]	4 (3,6)	3 (2,4)	0.34611
(70,80]	12 (-8,31)	2 (-17,3)	0.35342
vig_bout
(18,30]	1 (0,1)	1 (0,1)	0.02386
(30,40]	0 (0,1)	0 (0,1)	<0.001
(40,50]	0 (0,0)	1 (1,1)	0.19637
(50,60]	0 (0,0)	0 (0,1)	0.07122
(60,70]	0 (0,0)	0 (0,0)	0.34611
(70,80]	5 (-5,16)	0 (-10,0)	0.35342
x0_5
(18,30]	916 (902,931)	877 (863,895)	0.89457
(30,40]	874 (863,885)	858 (847,878)	0.52443
(40,50]	870 (849,890)	863 (843,882)	0.08212
(50,60]	875 (866,883)	900 (891,917)	0.26054
(60,70]	910 (893,926)	953 (937,972)	0.58360
(70,80]	956 (934,977)	981 (960,1005)	0.32919
x5_10
(18,30]	27 (26,28)	26 (25,27)	0.89457
(30,40]	28 (26,30)	25 (24,26)	0.52443
(40,50]	30 (29,32)	28 (27,30)	0.08212
(50,60]	31 (29,32)	28 (26,30)	0.26054
(60,70]	31 (30,32)	30 (28,31)	0.58360
(70,80]	32 (30,34)	29 (27,30)	0.32919
x10_25
(18,30]	50 (48,51)	48 (46,50)	<0.001
(30,40]	51 (49,54)	47 (45,49)	0.24128
(40,50]	54 (52,57)	52 (49,54)	0.67543
(50,60]	55 (53,58)	49 (46,52)	0.01238
(60,70]	54 (52,56)	51 (49,53)	<0.001
(70,80]	55 (51,58)	50 (47,53)	0.05672
x25_50
(18,30]	51 (49,52)	49 (48,51)	<0.001
(30,40]	52 (50,55)	50 (47,51)	0.24128
(40,50]	55 (52,58)	53 (50,55)	0.67543
(50,60]	56 (54,58)	49 (47,53)	0.01238
(60,70]	54 (52,55)	49 (47,50)	<0.001
(70,80]	53 (50,56)	50 (47,52)	0.05672
x50_75
(18,30]	36 (34,37)	35 (34,36)	0.11573
(30,40]	37 (36,39)	35 (34,37)	0.02764
(40,50]	38 (36,40)	37 (35,39)	0.06691
(50,60]	39 (38,41)	35 (34,37)	0.04701
(60,70]	38 (36,39)	33 (32,34)	0.07384
(70,80]	36 (34,38)	34 (33,36)	0.02859
x75_100
(18,30]	28 (27,29)	28 (27,29)	0.11573
(30,40]	30 (28,31)	28 (27,29)	0.02764
(40,50]	30 (28,32)	29 (28,30)	0.06691
(50,60]	31 (30,32)	28 (27,29)	0.04701
(60,70]	30 (29,30)	26 (25,27)	0.07384
(70,80]	29 (27,30)	27 (25,28)	0.02859
x100_150
(18,30]	43 (41,44)	43 (42,45)	0.15436
(30,40]	46 (44,48)	43 (41,44)	0.03387
(40,50]	47 (44,49)	45 (43,47)	0.12130
(50,60]	48 (47,49)	42 (41,44)	0.01598
(60,70]	46 (45,48)	40 (38,41)	0.00808
(70,80]	44 (42,47)	40 (38,42)	0.03169
x150_200
(18,30]	32 (31,33)	33 (32,34)	0.15436
(30,40]	35 (34,36)	32 (31,34)	0.03387
(40,50]	36 (34,37)	33 (32,35)	0.12130
(50,60]	37 (36,38)	33 (32,34)	0.01598
(60,70]	36 (34,37)	30 (29,31)	0.00808
(70,80]	34 (32,36)	30 (28,32)	0.03169
x200_500
(18,30]	109 (104,113)	112 (108,118)	0.18842
(30,40]	122 (119,125)	113 (110,118)	0.10493
(40,50]	124 (119,129)	113 (108,118)	0.28985
(50,60]	127 (122,131)	112 (108,116)	0.00586
(60,70]	123 (118,127)	104 (99,108)	<0.001
(70,80]	109 (102,116)	101 (93,106)	0.03914
x500_1e_03
(18,30]	77 (73,82)	85 (81,92)	0.18842
(30,40]	90 (85,95)	92 (87,98)	0.10493
(40,50]	87 (82,92)	87 (82,91)	0.28985
(50,60]	84 (80,87)	83 (79,87)	0.00586
(60,70]	74 (68,79)	70 (65,75)	<0.001
(70,80]	56 (51,60)	60 (55,66)	0.03914
x1e_03_2e_03
(18,30]	49 (46,52)	65 (62,71)	0.30178
(30,40]	55 (52,59)	76 (73,83)	0.13340
(40,50]	50 (47,54)	66 (63,72)	0.46801
(50,60]	43 (39,46)	58 (54,61)	<0.001
(60,70]	34 (29,39)	40 (35,45)	<0.001
(70,80]	20 (18,22)	29 (27,33)	0.08577
x2e_03_4e_03
(18,30]	19 (17,21)	30 (27,32)	0.30178
(30,40]	17 (15,19)	34 (32,38)	0.13340
(40,50]	16 (14,17)	27 (25,31)	0.46801
(50,60]	13 (11,15)	20 (18,21)	<0.001
(60,70]	10 (8,12)	13 (11,14)	<0.001
(70,80]	6 (4,7)	8 (7,9)	0.08577
x4e_03_6e_03
(18,30]	3 (2,4)	6 (5,7)	0.94802
(30,40]	2 (1,3)	5 (4,6)	0.02987
(40,50]	2 (2,3)	4 (3,5)	0.42978
(50,60]	2 (1,3)	4 (3,5)	<0.001
(60,70]	1 (0,2)	2 (1,2)	<0.001
(70,80]	5 (-4,14)	1 (-8,2)	0.03305
x6e_03_inf
(18,30]	1 (1,2)	2 (1,2)	0.94802
(30,40]	1 (0,1)	1 (1,2)	0.02987
(40,50]	0 (0,1)	1 (1,2)	0.42978
(50,60]	0 (0,1)	1 (1,1)	<0.001
(60,70]	0 (0,0)	0 (0,0)	<0.001
(70,80]	5 (-5,16)	0 (-10,0)	0.03305
total
(18,30]	254386 (240503,268269)	332115 (318232,350459)	0.49035
(30,40]	262405 (249548,275261)	356634 (343778,378669)	0.01703
(40,50]	247547 (234130,260964)	313775 (300358,334544)	0.42035
(50,60]	224743 (211646,237840)	270976 (257879,286446)	<0.001
(60,70]	190957 (177227,204686)	197972 (184243,212781)	<0.001
(70,80]	210066 (73864,346269)	157540 (21337,172983)	0.01666
mean
(18,30]	177 (167,186)	231 (221,243)	0.49035
(30,40]	182 (173,191)	248 (239,263)	0.01703
(40,50]	172 (163,181)	218 (209,232)	0.42035
(50,60]	156 (147,165)	188 (179,199)	<0.001
(60,70]	133 (123,142)	137 (128,148)	<0.001
(70,80]	146 (51,240)	109 (15,120)	0.01666
mvpa
(18,30]	23 (20,27)	37 (33,41)	0.57917
(30,40]	20 (17,22)	40 (37,44)	0.01839
(40,50]	18 (16,20)	32 (30,36)	0.09202
(50,60]	15 (12,17)	23 (21,26)	<0.001
(60,70]	11 (9,13)	14 (12,16)	<0.001
(70,80]	16 (-4,35)	9 (-11,10)	0.02074
mod
(18,30]	22 (19,25)	35 (32,38)	0.57917
(30,40]	19 (16,21)	39 (36,43)	0.01839
(40,50]	18 (15,20)	31 (29,34)	0.09202
(50,60]	14 (12,16)	23 (20,25)	<0.001
(60,70]	11 (9,13)	14 (12,16)	<0.001
(70,80]	10 (2,18)	9 (0,10)	0.02074
vig
(18,30]	1 (1,2)	2 (1,2)	0.16255
(30,40]	1 (0,1)	1 (1,2)	0.02064
(40,50]	0 (0,1)	1 (1,2)	0.01938
(50,60]	0 (0,1)	1 (1,1)	<0.001
(60,70]	0 (0,0)	0 (0,0)	<0.001
(70,80]	6 (-5,17)	0 (-11,0)	0.04391

2005-2006

Code

survey_design =
  df_small %>%
  filter(data_release_cycle == 4) %>%
  mutate(age_cat = cut(age_in_years_at_screening,
                       breaks=c(18, 30, seq(40, 80, 10)))) %>%
  mutate(weight = full_sample_2_year_mec_exam_weight,
         weight_norm = weight / mean(weight)) %>%
  ungroup()  %>%
  survey::svydesign(
    id = ~ masked_variance_pseudo_psu,
    strata = ~ masked_variance_pseudo_stratum,
    weights = ~ weight_norm,
    data = .,
    nest = TRUE
  )

calc_by_age =
  function(var, df) {
    formula = as.formula(paste("~", var))
    total_by_age_gender = svyby(formula,
                                ~ age_cat + gender,
                                df,
                                svymean) %>%
      rename(mean = contains(var)) %>%
      mutate(metric = var)
  }

means_df =
  map_dfr(.x = pa_old %>% select(contains("bout"), starts_with("x"), total, mean, mvpa, mod, vig) %>% colnames(),
          .f = calc_by_age, df = survey_design)


small = 
  means_df %>%
  pivot_wider(names_from = gender, values_from = c(mean, se)) %>% 
  mutate(val_female = paste0(round(mean_Female), " (", round(mean_Female - 1.96 * se_Female), ",", round(mean_Female + 1.96 * se_Female), ")"),
         val_male = paste0(round(mean_Male), " (", round(mean_Male - 1.96 * se_Female), ",", round(mean_Male + 1.96 * se_Male), ")"))

# test = svyttest(mvpa_montoye ~ gender, survey_design)
pvals = c()
for(i in 1:nrow(small)){
  row = means_df[i,]
  cat = row$age_cat
  var = row$metric
  survey_design =
    df_small %>%
    filter(age_in_years_at_screening >= 18 & data_release_cycle == 4) %>%
    mutate(age_cat = cut(age_in_years_at_screening,
                         breaks=c(18, 30, seq(40, 80, 10)))) %>%
    filter(age_cat == cat) %>%
    mutate(weight = full_sample_2_year_mec_exam_weight / 2,
           weight_norm = weight / mean(weight)) %>%
    ungroup()  %>%
    survey::svydesign(
      id = ~ masked_variance_pseudo_psu,
      strata = ~ masked_variance_pseudo_stratum,
      weights = ~ weight_norm,
      data = .,
      nest = TRUE
    )
  formula = as.formula(paste(var,  "~ gender"))
  test = svyttest(formula, survey_design)
  pvals = append(pvals, unname(test$p.value))
}

small %>% 
  select(metric, age_cat, Female = val_female, Male = val_male, mean_Female, mean_Male) %>% 
  bind_cols(p = pvals) %>% 
  mutate(pval = if_else(p < 0.001, "<0.001", format.pval(p, digits = 2))) %>% 
  # mutate(metric = factor(metric, 
  #                        labels = c("MVPA (min)", 
  #                                   "VPA (min)", "Time <1000",
  #                                   "Time [1000,2000]", "Time (2000,3000]", "Time (3000,4000]", "Time (4000,5000]", "Time (5000,6000]", "Time (6000,7000]", "Time (7000,8000]", "Time (8000,9000]", "Time >9000", "Total AC", "Mean AC", "Total MIMS"))) %>% 
  group_by(metric) %>% 
  gt() %>% 
  cols_hide(c("mean_Female", "mean_Male", "p")) %>% 
  tab_style(
    # Style to apply when conditions are met
    style = cell_fill(color = "lightpink"),
    locations = cells_body(
      rows = mean_Female > mean_Male & p < 0.05
    )
  ) %>% 
  tab_style(
    # Style to apply when conditions are met
    style = cell_fill(color = "lightblue"),
    locations = cells_body(
      rows = mean_Female < mean_Male & p < 0.05
    )
  )

age_cat	Female	Male	pval
mvpa_bout
(18,30]	3 (2,4)	6 (5,7)	0.00501
(30,40]	4 (3,5)	4 (3,5)	0.74898
(40,50]	4 (3,5)	4 (3,5)	0.82775
(50,60]	3 (2,4)	4 (3,5)	0.23840
(60,70]	2 (2,3)	4 (4,6)	0.07196
(70,80]	3 (2,4)	2 (1,3)	0.38992
vig_bout
(18,30]	0 (0,1)	0 (0,1)	0.00501
(30,40]	0 (0,1)	1 (0,1)	0.74898
(40,50]	0 (0,1)	1 (1,1)	0.82775
(50,60]	0 (0,0)	0 (0,1)	0.23840
(60,70]	0 (0,0)	0 (0,0)	0.07196
(70,80]	0 (0,0)	0 (0,0)	0.38992
x0_5
(18,30]	911 (900,923)	877 (865,890)	0.52682
(30,40]	884 (874,895)	851 (841,870)	0.72103
(40,50]	873 (855,891)	877 (859,898)	0.27731
(50,60]	880 (863,898)	873 (856,899)	0.06383
(60,70]	921 (903,939)	917 (898,936)	0.24140
(70,80]	940 (912,969)	987 (958,999)	0.21847
x5_10
(18,30]	28 (27,29)	26 (25,27)	0.52682
(30,40]	29 (28,30)	28 (27,29)	0.72103
(40,50]	29 (28,31)	27 (26,29)	0.27731
(50,60]	30 (29,32)	30 (28,31)	0.06383
(60,70]	31 (29,32)	30 (28,31)	0.24140
(70,80]	32 (30,34)	31 (29,32)	0.21847
x10_25
(18,30]	51 (50,53)	48 (46,50)	0.00132
(30,40]	53 (51,55)	52 (50,53)	0.00376
(40,50]	53 (51,56)	50 (48,53)	0.78499
(50,60]	55 (53,57)	53 (51,56)	0.59360
(60,70]	53 (50,56)	53 (50,56)	0.61409
(70,80]	57 (54,59)	53 (50,56)	0.00473
x25_50
(18,30]	51 (50,52)	50 (49,51)	0.00132
(30,40]	54 (51,56)	53 (51,55)	0.00376
(40,50]	54 (52,56)	52 (50,54)	0.78499
(50,60]	54 (52,56)	54 (52,57)	0.59360
(60,70]	53 (50,56)	53 (50,56)	0.61409
(70,80]	56 (53,58)	52 (49,54)	0.00473
x50_75
(18,30]	36 (35,36)	36 (35,37)	0.01027
(30,40]	38 (37,39)	37 (36,39)	0.23841
(40,50]	38 (37,39)	36 (35,38)	0.06667
(50,60]	38 (37,39)	38 (37,39)	0.30456
(60,70]	37 (35,39)	36 (34,39)	0.27720
(70,80]	39 (36,41)	35 (33,37)	0.29542
x75_100
(18,30]	28 (27,28)	28 (27,29)	0.01027
(30,40]	30 (28,31)	29 (28,31)	0.23841
(40,50]	30 (29,31)	28 (27,30)	0.06667
(50,60]	30 (29,31)	30 (29,31)	0.30456
(60,70]	29 (27,30)	28 (26,30)	0.27720
(70,80]	30 (29,31)	27 (25,28)	0.29542
x100_150
(18,30]	43 (42,44)	42 (41,44)	0.00390
(30,40]	45 (44,47)	44 (43,47)	0.46062
(40,50]	45 (44,47)	43 (41,45)	0.08025
(50,60]	45 (44,47)	45 (43,47)	0.34405
(60,70]	44 (42,46)	43 (40,45)	0.69794
(70,80]	47 (44,49)	39 (37,41)	0.08751
x150_200
(18,30]	32 (31,33)	32 (31,33)	0.00390
(30,40]	34 (33,35)	33 (32,35)	0.46062
(40,50]	34 (33,36)	32 (30,34)	0.08025
(50,60]	35 (34,36)	33 (32,35)	0.34405
(60,70]	34 (33,36)	32 (31,34)	0.69794
(70,80]	35 (33,37)	29 (27,30)	0.08751
x200_500
(18,30]	110 (106,113)	107 (104,112)	0.25684
(30,40]	115 (111,119)	115 (111,120)	0.65398
(40,50]	118 (113,123)	109 (104,114)	0.20979
(50,60]	123 (118,127)	114 (109,119)	0.94003
(60,70]	119 (115,123)	109 (105,115)	0.86345
(70,80]	114 (104,123)	94 (84,98)	0.05127
x500_1e_03
(18,30]	81 (77,84)	84 (81,89)	0.25684
(30,40]	84 (79,88)	92 (87,97)	0.65398
(40,50]	88 (83,94)	86 (80,89)	0.20979
(50,60]	86 (81,91)	85 (80,91)	0.94003
(60,70]	75 (71,78)	76 (72,82)	0.86345
(70,80]	59 (53,66)	56 (49,61)	0.05127
x1e_03_2e_03
(18,30]	51 (48,54)	69 (66,75)	0.80340
(30,40]	53 (49,57)	72 (68,76)	0.65742
(40,50]	55 (52,57)	66 (63,70)	0.19569
(50,60]	48 (44,52)	60 (56,65)	0.89445
(60,70]	35 (32,38)	47 (43,50)	0.52245
(70,80]	24 (20,28)	28 (24,32)	0.03411
x2e_03_4e_03
(18,30]	16 (15,17)	33 (32,37)	0.80340
(30,40]	17 (15,18)	29 (27,31)	0.65742
(40,50]	18 (17,19)	28 (26,31)	0.19569
(50,60]	13 (11,15)	21 (19,24)	0.89445
(60,70]	9 (7,10)	15 (13,16)	0.52245
(70,80]	7 (5,8)	7 (6,9)	0.03411
x4e_03_6e_03
(18,30]	2 (2,3)	6 (6,8)	0.45105
(30,40]	3 (2,4)	4 (3,5)	0.52543
(40,50]	3 (2,3)	4 (3,4)	0.23567
(50,60]	2 (1,2)	3 (2,4)	0.58007
(60,70]	1 (1,1)	3 (2,4)	0.41695
(70,80]	1 (0,1)	1 (0,2)	0.00244
x6e_03_inf
(18,30]	1 (1,1)	2 (2,2)	0.45105
(30,40]	1 (1,1)	1 (1,2)	0.52543
(40,50]	1 (1,1)	1 (1,2)	0.23567
(50,60]	0 (0,0)	1 (1,1)	0.58007
(60,70]	0 (0,0)	0 (0,1)	0.41695
(70,80]	0 (0,0)	0 (0,0)	0.00244
total
(18,30]	244468 (233741,255194)	345691 (334965,367519)	0.59801
(30,40]	256310 (246017,266603)	331941 (321648,344872)	0.47110
(40,50]	264732 (254030,275434)	313103 (302401,335408)	0.16915
(50,60]	230002 (213559,246445)	276450 (260007,300438)	0.56800
(60,70]	186868 (174751,198986)	228355 (216237,240040)	0.16732
(70,80]	151700 (133850,169550)	152349 (134500,161934)	<0.001
mean
(18,30]	170 (162,177)	240 (233,255)	0.59801
(30,40]	178 (171,185)	231 (223,239)	0.47110
(40,50]	184 (176,191)	217 (210,233)	0.16915
(50,60]	160 (148,171)	192 (181,209)	0.56800
(60,70]	130 (121,138)	159 (150,167)	0.16732
(70,80]	105 (93,118)	106 (93,112)	<0.001
mvpa
(18,30]	19 (18,21)	40 (39,45)	0.94371
(30,40]	20 (18,22)	33 (32,36)	0.26065
(40,50]	21 (19,23)	32 (30,36)	0.04576
(50,60]	14 (12,17)	24 (22,28)	0.11235
(60,70]	10 (8,12)	17 (16,20)	0.02303
(70,80]	7 (5,9)	9 (7,11)	<0.001
mod
(18,30]	18 (17,20)	38 (37,43)	0.94371
(30,40]	19 (18,21)	32 (30,34)	0.26065
(40,50]	20 (18,22)	31 (29,34)	0.04576
(50,60]	14 (12,17)	23 (21,27)	0.11235
(60,70]	10 (8,11)	17 (15,19)	0.02303
(70,80]	7 (5,9)	8 (6,10)	<0.001
vig
(18,30]	1 (1,1)	2 (2,2)	0.33546
(30,40]	1 (1,1)	1 (1,2)	0.90075
(40,50]	1 (1,1)	1 (1,2)	0.02948
(50,60]	0 (0,0)	1 (1,1)	0.02907
(60,70]	0 (0,0)	0 (0,1)	0.01403
(70,80]	0 (0,0)	0 (0,0)	0.00142

Plots

Distributions

Code

pa_old %>% 
  left_join(demo) %>% 
  filter(age_in_years_at_screening >= 18) %>% 
  pivot_longer(cols = c(mvpa, mvpa_bout, vig, vig_bout, mean)) %>% 
  mutate(name = factor(name, labels = c("Mean AC", "MVPA (min)", "MVPA (10 min bout min)", "VPA (min)", "VPA (10 min bout min)"))) %>% 
  ggplot(aes(x = value, color = gender)) + 
  geom_density() + 
  facet_wrap(.~name, scales = "free") + 
  scale_color_manual(values = c(col1, col2), name = "")+
  theme(legend.position = c(0.9, 0.2)) + 
  labs(x = "Value", y = "Density")

We see that the data are very highly skewed, so we can Winsorize at the 99th percentile:

Code

pa_old_win = 
  pa_old %>% 
    left_join(demo) %>% 
  filter(age_in_years_at_screening >= 18) %>% 
  mutate(across(c(mvpa, mod, vig, mean, total, mvpa_bout, vig_bout),
                ~DescTools::Winsorize(.x, quantile(.x, probs = c(0, 0.99)))))

pa_old_win %>% 
  pivot_longer(cols = c(mvpa, mvpa_bout, vig, vig_bout, mean)) %>% 
  mutate(name = factor(name, labels = c("Mean AC", "MVPA (min)", "MVPA (10 min bout min)", "VPA (min)", "VPA (10 min bout min)"))) %>% 
  ggplot(aes(x = value, color = gender)) + 
  geom_density() + 
  facet_wrap(.~name, scales = "free") + 
  scale_color_manual(values = c(col1, col2), name = "")+
  theme(legend.position = c(0.9, 0.2)) + 
  labs(x = "Value", y = "Density")

Smoothed, unweighted

Code

pa_old_win %>%
  pivot_longer(cols = c(contains("bout"), vig, mvpa)) %>% 
  mutate(name = factor(name, labels = c("MVPA", "MVPA 10 min bout", "VPA", "VPA 10 min bout"))) %>%
  ggplot(aes(x = age_in_years_at_screening, y = value, color = gender)) +
  geom_jitter(alpha = .1, size = .7, width = .25) +
  geom_smooth(se = F) +
  facet_wrap(name~., scales = "free_y") + 
  labs(x = "Age", y = "Minutes", color = "")+
  scale_color_manual(values = c(col1, col2)) +
  theme(legend.position = c(0.9, 0.9),
        legend.title = element_blank())

Code

# pa_old_win %>%
#   pivot_longer(cols = c("mvpa", "vig", "mod")) %>% 
#   mutate(name = factor(name, labels = c("Mod", "MVPA", "VPA"))) %>% 
#   ggplot(aes(x = age_in_years_at_screening, y = value, color = gender)) +
#   geom_jitter(alpha = .1, size = .7, width = .25) +
#   geom_smooth(se = F) +
#   facet_grid(name~data_release_cycle, scales = "free_y") + 
#   labs(x = "Age", y = "Minutes", color = "")+
#   scale_color_manual(values = c(col1, col2)) +
#   theme(legend.position = c(0.9, 0.9),
#         legend.title = element_blank())

pa_old_win %>%
  pivot_longer(cols = starts_with("x")) %>% 
  mutate(name = factor(
    name,
    levels = c(
      "x0_5",
      "x5_10",
      "x10_25",
      "x25_50",
      "x50_75",
      "x75_100" ,
      "x100_150",
      "x150_200",
      "x200_500" ,
      "x500_1e_03",
      "x1e_03_2e_03",
      "x2e_03_4e_03",
      "x4e_03_6e_03",
      "x6e_03_inf"
    )
  )) %>% 
  ggplot(aes(x = age_in_years_at_screening, y = value, color = gender)) +
  geom_jitter(alpha = .1, size = .7, width = .25) +
  geom_smooth(se = F) +
  facet_wrap(name~., scales = "free_y") + 
  labs(x = "Age", y = "Minutes", color = "")+
  scale_color_manual(values = c(col1, col2)) +
  theme(legend.position = "bottom",
        legend.title = element_blank())

Smoothed, weighted

Code

survey_design =
  df_small %>%
  mutate(across(c(mvpa, mod, vig, mean, total, starts_with("x"), mvpa_bout, vig_bout),
                ~DescTools::Winsorize(.x, quantile(.x, probs = c(0, 0.99))))) %>% 
  filter(age_in_years_at_screening >= 18) %>%
  mutate(weight = full_sample_2_year_mec_exam_weight / 2,
         weight_norm = weight / mean(weight)) %>%
  ungroup()  %>%
  survey::svydesign(
    id = ~ masked_variance_pseudo_psu,
    strata = ~ masked_variance_pseudo_stratum,
    weights = ~ weight_norm,
    data = .,
    nest = TRUE
  )


calc_by_age =
  function(var, df) {
    # var = "total_oaksteps"
    formula = as.formula(paste("~", var))
    total_by_age_gender = svyby(formula,
                                ~ age_in_years_at_screening + gender,
                                df,
                                svymean) %>%
      rename(mean = contains(var)) %>%
      mutate(metric = var)
  }

means_df =
  map_dfr(.x = pa_old %>% select(contains("bout"), mvpa, mod, vig, mean, total) %>% colnames(),
          .f = calc_by_age, df = survey_design)

# get smooths for means and confidence intervals
models = means_df %>%
  mutate(lb = mean - 1.96 * se,
         ub = mean + 1.96 * se) %>%
  tidyr::nest(data = -c(metric, gender)) %>%
  dplyr::mutate(
    # Perform loess calculation on each group
    m = purrr::map(data, loess,
                   formula = mean ~ age_in_years_at_screening, span = .75),
    # Retrieve the fitted values from each model
    fitted_mean = purrr::map(m, `[[`, "fitted"),
    l = purrr::map(data, loess,
                   formula = lb ~ age_in_years_at_screening, span = .75),
    # Retrieve the fitted values from each model
    fitted_lb = purrr::map(l, `[[`, "fitted"),
    u = purrr::map(data, loess,
                   formula = ub ~ age_in_years_at_screening, span = .75),
    # Retrieve the fitted values from each model
    fitted_ub = purrr::map(u, `[[`, "fitted")
  )

# Apply fitted y's as a new column
results = models %>%
  dplyr::select(-m, -l, -u) %>%
  tidyr::unnest(cols = c(data, contains("fitted")))

results %>%
  mutate(metric = factor(metric, labels = c("Mean AC", "MPA (min)", "MVPA (min)", "MVPA 10 min bout (min)", "Total AC", "VPA (min)", "VPA 10 min bout (min)"))) %>%
    ggplot(aes(x = age_in_years_at_screening, y = fitted_mean,
             ymin = fitted_lb, ymax = fitted_ub, color = gender, fill = gender)) +
  facet_wrap(. ~ metric, scales = "free_y") +
  geom_line(linewidth = 1) +
  geom_ribbon(alpha = .2, linetype = 0) +
  theme_light() +
  theme(legend.position = "bottom",
        legend.title = element_blank(),
        text =element_text(size = 14),
        strip.text = element_text(size = 12),
        panel.grid.minor = element_blank()) +
  labs(x = "Age (years)", y = "Mean Daily Minutes",
       title = "Smoothed survey weighted mean daily minutes by age ")+
  scale_x_continuous(breaks=seq(20,80,10)) +
  scale_y_continuous(breaks=seq(0,600,30)) + 
  scale_color_manual(values = c(col1, col2))+
  scale_fill_manual(values = c(col1, col2))

Code

means_df =
  map_dfr(.x = pa_old %>% select(starts_with("x")) %>% colnames(),
          .f = calc_by_age, df = survey_design)

# get smooths for means and confidence intervals
models = means_df %>%
  mutate(lb = mean - 1.96 * se,
         ub = mean + 1.96 * se) %>%
  tidyr::nest(data = -c(metric, gender)) %>%
  dplyr::mutate(
    # Perform loess calculation on each group
    m = purrr::map(data, loess,
                   formula = mean ~ age_in_years_at_screening, span = .75),
    # Retrieve the fitted values from each model
    fitted_mean = purrr::map(m, `[[`, "fitted"),
    l = purrr::map(data, loess,
                   formula = lb ~ age_in_years_at_screening, span = .75),
    # Retrieve the fitted values from each model
    fitted_lb = purrr::map(l, `[[`, "fitted"),
    u = purrr::map(data, loess,
                   formula = ub ~ age_in_years_at_screening, span = .75),
    # Retrieve the fitted values from each model
    fitted_ub = purrr::map(u, `[[`, "fitted")
  )

# Apply fitted y's as a new column
results = models %>%
  dplyr::select(-m, -l, -u) %>%
  tidyr::unnest(cols = c(data, contains("fitted")))

results %>%
   mutate(metric = factor(
    metric,
    levels = c(
      "x0_5",
      "x5_10",
      "x10_25",
      "x25_50",
      "x50_75",
      "x75_100" ,
      "x100_150",
      "x150_200",
      "x200_500" ,
      "x500_1e_03",
      "x1e_03_2e_03",
      "x2e_03_4e_03",
      "x4e_03_6e_03",
      "x6e_03_inf"
    )
  )) %>% 
    ggplot(aes(x = age_in_years_at_screening, y = fitted_mean,
             ymin = fitted_lb, ymax = fitted_ub, color = gender, fill = gender)) +
  facet_wrap(. ~ metric, scales = "free_y") +
  geom_line(linewidth = 1) +
  geom_ribbon(alpha = .2, linetype = 0) +
  theme_light() +
  theme(legend.position = "bottom",
        legend.title = element_blank(),
        text =element_text(size = 14),
        strip.text = element_text(size = 12),
        panel.grid.minor = element_blank()) +
  labs(x = "Age (years)", y = "Mean Daily Minutes",
       title = "Smoothed survey weighted mean daily minutes by age")+
  scale_x_continuous(breaks=seq(20,80,10)) +
  scale_color_manual(values = c(col1, col2))+
  scale_fill_manual(values = c(col1, col2))

Investigating outliers

Code

outliers = 
  pa_old %>% 
  left_join(demo) %>% 
  filter(age_in_years_at_screening >= 18) %>% 
  arrange(desc(total)) %>% 
  slice(1:10) 

# pa_old %>% 
#   arrange(desc(mvpa)) %>% 
#   slice(1:10) %>% 
#   select(SEQN, total, mean, mvpa)

# by day 
# pa_byday = readRDS(here::here("data", "accelerometry", 
#                               "summarized", "pa_summary_AC_C_D.rds")) %>% 
#   mutate(SEQN = as.character(SEQN)) 
# 
# pa_byday %>% 
#   filter(SEQN == "26729") %>% 
#   select(SEQN, PAXDAY, mod, vig, mvpa)

minute_level = readRDS(here::here("data", "accelerometry", "minute_level", "NHANES_1440_C_PAXINTEN.rds")) %>% 
  bind_rows(readRDS(here::here("data", "accelerometry", "minute_level", "NHANES_1440_D_PAXINTEN.rds")))

for(sub in outliers$SEQN){
  temp = 
    minute_level %>% 
    filter(SEQN == sub)
  
  mvpa = round(outliers %>% filter(SEQN == sub) %>% pull(mvpa), 0)
  age = demo %>% filter(SEQN == sub) %>% pull(age_in_years_at_screening)
  sex = demo %>% filter(SEQN == sub) %>% pull(gender)
  p = 
    temp %>% 
    pivot_longer(cols = starts_with("min")) %>% 
    mutate(m = sub(".*min\\_", "", name) %>% as.numeric()) %>%
    ggplot(aes(x = m, y = value)) + 
    geom_line() + 
    facet_grid(PAXDAY~.) + 
    geom_hline(aes(yintercept = 2020), color = "lightblue", linetype = "dashed") + 
    scale_x_continuous(breaks=seq(0,1440,120), labels =seq(0,24,2)) + 
    labs(x = "Hour of Day", y = "Counts", title = paste0("SEQN = ", sub, " MVPA = ", mvpa),
         subtitle = paste0("Age = ", age, " Sex = ", sex))
  
  print(p)
}

NHANES 2011-2014

Code

# load packages 


pa = readRDS(here::here("data", "accelerometry", 
                              "summarized", "pa_summary_AC.rds")) %>% 
  mutate(SEQN = as.character(SEQN)) %>% 
  group_by(SEQN) %>%
  summarize(across(-PAXDAYM, ~mean(.x, na.rm = TRUE))) 

days = read_csv(here::here("data", "accelerometry", "inclusion_summary.csv.gz"))

wear_gh = 
  days %>% 
  group_by(SEQN) %>% 
  summarize(wear_mins = mean(non_flag_wear)) %>% 
  mutate(SEQN = as.character(SEQN))

demo_gh = readRDS(here::here("data", "demographics",  "covariates_accel_mortality_df.rds")) %>% 
  select(-contains("peak"), -contains("total_"))

Survey Weighted Physical Activity Summaries by Sex

Code

df_small_gh = pa %>% 
  left_join(demo_gh) %>% 
  left_join(wear_gh) %>% 
  filter(age_in_years_at_screening >= 18) %>%
  select(SEQN,
         data_release_cycle,
         gender,
         wear_mins,
         masked_variance_pseudo_psu,
         masked_variance_pseudo_stratum,
         full_sample_2_year_interview_weight,
         full_sample_2_year_mec_exam_weight,
         age_in_years_at_screening,
         total, 
         mean,
         mvpa_montoye, 
         mod_montoye,
         vig_montoye,
         mvpa_montoye_bout,
         vig_montoye_bout,
         starts_with("x"),
         missing_ac = na) 
labs = c("SEQN", "Data release cycle", "Gender",  "Wear minutes",
         "Pseudo PSU",
         "Psueudo stratum",
         "2 yr int weight", "2 yr exam weight",
         "Age", "Total AC", "Mean AC", "MVPA minutes",  "Moderate minutes", "Vigorous minutes",
         "MVPA (10 min bout)",
         "Vig (10 min bout)",
         "Time <100", "Time [100,500]","Time (500,1000]","Time (1000,2000]", "Time (2000,3000]", "Time (3000,4000]", "Time (4000,5000]", "Time (5000,6000]", "Time (6000,7000]",  "Time (7000,8000]", "Time (8000,9000]", "Time >9000", "Missing AC") 

names(labs) = colnames(df_small_gh)

df_small_gh =
  df_small_gh %>%
  labelled::set_variable_labels(!!!labs)

df_svy =
  df_small_gh %>% 
  mutate(WTMEC4YR = full_sample_2_year_mec_exam_weight/2,
         WTMEC4YR_norm = WTMEC4YR/mean(WTMEC4YR, na.rm = TRUE)) %>%
  svydesign(ids = ~masked_variance_pseudo_psu, weights = ~WTMEC4YR_norm,
            strata = ~masked_variance_pseudo_stratum, data=., nest=TRUE)

# survey weighted table
df_svy %>%
  tbl_svysummary(
    by = gender,
    include = -c(SEQN, data_release_cycle, masked_variance_pseudo_psu, masked_variance_pseudo_stratum, full_sample_2_year_interview_weight, full_sample_2_year_mec_exam_weight, WTMEC4YR, WTMEC4YR_norm),
    statistic = list(
      all_continuous() ~ "{mean} ({sd})",
      all_categorical() ~ "{n} ({p}%)"
    ),
    digits = all_continuous() ~ 0,
    missing_text = "Missing",
  ) %>%
  add_overall() %>%
  add_p()

Characteristic	Overall N = 9,298 1	Female N = 4,876 1	Male N = 4,422 1	p-value 2
Wear minutes	1,176 (137)	1,177 (135)	1,174 (140)	0.4
Age	47 (17)	48 (18)	47 (17)	<0.001
Total AC	2,557,250 (825,267)	2,648,375 (817,869)	2,456,752 (821,806)	<0.001
Mean AC	1,778 (574)	1,841 (569)	1,708 (572)	<0.001
MVPA minutes	299 (120)	314 (118)	282 (119)	<0.001
Moderate minutes	122 (41)	125 (40)	118 (43)	<0.001
Vigorous minutes	137 (84)	149 (86)	125 (81)	<0.001
MVPA (10 min bout)	99 (84)	108 (85)	89 (82)	<0.001
Vig (10 min bout)	32 (44)	35 (45)	28 (41)	<0.001
Time <100	574 (105)	580 (102)	568 (108)	<0.001
Time [100,500]	129 (39)	121 (35)	137 (40)	<0.001
Time (500,1000]	119 (29)	111 (26)	127 (29)	<0.001
Time (1000,2000]	160 (34)	155 (31)	166 (35)	<0.001
Time (2000,3000]	112 (25)	112 (24)	112 (27)	0.5
Time (3000,4000]	91 (24)	92 (23)	90 (25)	0.015
Time (4000,5000]	77 (25)	79 (24)	75 (26)	<0.001
Time (5000,6000]	62 (25)	65 (25)	59 (26)	<0.001
Time (6000,7000]	45 (23)	48 (23)	42 (23)	<0.001
Time (7000,8000]	30 (20)	32 (20)	27 (19)	<0.001
Time (8000,9000]	18 (14)	20 (15)	15 (14)	<0.001
Time >9000	23 (27)	25 (28)	20 (25)	<0.001
Missing AC	20 (21)	22 (22)	18 (20)	0.11
Missing	6,809	3,736	3,073
1 Mean (SD)
2 Design-based KruskalWallis test

Stratified by wave

2011-2012

Code

df_svy =
  df_small_gh %>% 
  filter(data_release_cycle == 7) %>%
  mutate(WTMEC4YR = full_sample_2_year_mec_exam_weight,
         WTMEC4YR_norm = WTMEC4YR/mean(WTMEC4YR, na.rm = TRUE)) %>%
  svydesign(ids = ~masked_variance_pseudo_psu, weights = ~WTMEC4YR_norm,
            strata = ~masked_variance_pseudo_stratum, data=., nest=TRUE)

# survey weighted table
df_svy %>%
  tbl_svysummary(
    by = gender,
    include = -c(SEQN, data_release_cycle, masked_variance_pseudo_psu, masked_variance_pseudo_stratum, full_sample_2_year_interview_weight, full_sample_2_year_mec_exam_weight, WTMEC4YR, WTMEC4YR_norm),
    statistic = list(
      all_continuous() ~ "{mean} ({sd})",
      all_categorical() ~ "{n} ({p}%)"
    ),
    digits = all_continuous() ~ 0,
    missing_text = "Missing",
  ) %>%
  add_overall() %>%
  add_p()

Characteristic	Overall N = 4,609 1	Female N = 2,390 1	Male N = 2,219 1	p-value 2
Wear minutes	1,178 (136)	1,180 (131)	1,176 (141)	0.3
Age	47 (17)	48 (17)	46 (17)	<0.001
Total AC	2,565,566 (818,326)	2,655,689 (812,943)	2,468,463 (813,162)	<0.001
Mean AC	1,784 (569)	1,846 (565)	1,716 (566)	<0.001
MVPA minutes	299 (119)	315 (119)	282 (117)	<0.001
Moderate minutes	122 (41)	126 (41)	118 (42)	<0.001
Vigorous minutes	137 (83)	149 (85)	125 (79)	<0.001
MVPA (10 min bout)	99 (83)	108 (85)	89 (79)	<0.001
Vig (10 min bout)	32 (43)	35 (45)	27 (40)	<0.001
Time <100	570 (104)	576 (101)	564 (107)	0.013
Time [100,500]	129 (39)	121 (36)	137 (40)	<0.001
Time (500,1000]	119 (28)	111 (26)	128 (28)	<0.001
Time (1000,2000]	162 (34)	156 (32)	168 (35)	<0.001
Time (2000,3000]	113 (26)	112 (24)	113 (27)	0.7
Time (3000,4000]	91 (24)	93 (23)	90 (25)	0.071
Time (4000,5000]	77 (25)	80 (25)	75 (26)	0.003
Time (5000,6000]	62 (25)	65 (25)	59 (25)	<0.001
Time (6000,7000]	45 (23)	48 (24)	42 (23)	<0.001
Time (7000,8000]	30 (19)	32 (20)	27 (18)	<0.001
Time (8000,9000]	18 (14)	20 (15)	15 (13)	<0.001
Time >9000	23 (26)	25 (26)	20 (26)	<0.001
Missing AC	19 (21)	21 (21)	18 (20)	0.3
Missing	3,404	1,855	1,549
1 Mean (SD)
2 Design-based KruskalWallis test

2013-2014

Code

df_svy =
  df_small_gh %>% 
  filter(data_release_cycle == 8) %>%
  mutate(WTMEC4YR = full_sample_2_year_mec_exam_weight,
         WTMEC4YR_norm = WTMEC4YR/mean(WTMEC4YR, na.rm = TRUE)) %>%
  svydesign(ids = ~masked_variance_pseudo_psu, weights = ~WTMEC4YR_norm,
            strata = ~masked_variance_pseudo_stratum, data=., nest=TRUE)

# survey weighted table
df_svy %>%
  tbl_svysummary(
    by = gender,
    include = -c(SEQN, data_release_cycle, masked_variance_pseudo_psu, masked_variance_pseudo_stratum, full_sample_2_year_interview_weight, full_sample_2_year_mec_exam_weight, WTMEC4YR, WTMEC4YR_norm),
    statistic = list(
      all_continuous() ~ "{mean} ({sd})",
      all_categorical() ~ "{n} ({p}%)"
    ),
    digits = all_continuous() ~ 0,
    missing_text = "Missing",
  ) %>%
  add_overall() %>%
  add_p()

Characteristic	Overall N = 4,689 1	Female N = 2,487 1	Male N = 2,202 1	p-value 2
Wear minutes	1,173 (139)	1,173 (139)	1,173 (139)	>0.9
Age	47 (18)	48 (18)	47 (17)	0.015
Total AC	2,548,867 (832,210)	2,641,164 (822,799)	2,444,654 (830,654)	<0.001
Mean AC	1,772 (579)	1,836 (572)	1,700 (578)	<0.001
MVPA minutes	298 (121)	312 (118)	282 (122)	<0.001
Moderate minutes	121 (41)	124 (39)	118 (43)	<0.001
Vigorous minutes	137 (85)	149 (86)	125 (83)	<0.001
MVPA (10 min bout)	99 (85)	108 (85)	90 (85)	<0.001
Vig (10 min bout)	32 (45)	36 (46)	28 (43)	<0.001
Time <100	578 (106)	584 (103)	572 (108)	0.002
Time [100,500]	128 (38)	120 (34)	138 (40)	<0.001
Time (500,1000]	118 (29)	110 (26)	126 (30)	<0.001
Time (1000,2000]	159 (33)	154 (31)	164 (34)	<0.001
Time (2000,3000]	111 (25)	111 (24)	112 (26)	0.5
Time (3000,4000]	90 (24)	91 (23)	90 (25)	0.082
Time (4000,5000]	77 (25)	78 (24)	75 (26)	0.001
Time (5000,6000]	62 (25)	64 (24)	59 (26)	<0.001
Time (6000,7000]	45 (23)	48 (23)	42 (24)	<0.001
Time (7000,8000]	30 (20)	32 (20)	27 (19)	<0.001
Time (8000,9000]	18 (15)	20 (15)	15 (14)	<0.001
Time >9000	23 (27)	26 (29)	19 (25)	<0.001
Missing AC	20 (21)	22 (23)	18 (20)	0.2
Missing	3,404	1,881	1,523
1 Mean (SD)
2 Design-based KruskalWallis test

Survey-Weighted Physical Activity Summaries by Sex and Age

Code

survey_design =
  df_small_gh %>%
  mutate(age_cat = cut(age_in_years_at_screening,
                       breaks=c(18, 30, seq(40, 80, 10)))) %>%
  mutate(weight = full_sample_2_year_mec_exam_weight / 2,
         weight_norm = weight / mean(weight)) %>%
  ungroup()  %>%
  survey::svydesign(
    id = ~ masked_variance_pseudo_psu,
    strata = ~ masked_variance_pseudo_stratum,
    weights = ~ weight_norm,
    data = .,
    nest = TRUE
  )

calc_by_age =
  function(var, df) {
    formula = as.formula(paste("~", var))
    total_by_age_gender = svyby(formula,
                                ~ age_cat + gender,
                                df,
                                svymean) %>%
      rename(mean = contains(var)) %>%
      mutate(metric = var)
  }

means_df =
  map_dfr(.x = pa %>% select(total, 
         mean,
         mvpa_montoye, 
         mod_montoye,
         vig_montoye,
         mvpa_montoye_bout,
         vig_montoye_bout,
         starts_with("x")) %>% colnames(),
          .f = calc_by_age, df = survey_design)


small = 
  means_df %>%
  pivot_wider(names_from = gender, values_from = c(mean, se)) %>% 
  mutate(val_female = paste0(round(mean_Female), " (", round(mean_Female - 1.96 * se_Female), ",", round(mean_Female + 1.96 * se_Female), ")"),
         val_male = paste0(round(mean_Male), " (", round(mean_Male - 1.96 * se_Female), ",", round(mean_Male + 1.96 * se_Male), ")"))

# test = svyttest(mvpa_montoye ~ gender, survey_design)
pvals = c()
for(i in 1:nrow(small)){
  row = means_df[i,]
  cat = row$age_cat
  var = row$metric
  survey_design =
    df_small_gh %>%
    mutate(age_cat = cut(age_in_years_at_screening,
                         breaks=c(18, 30, seq(40, 80, 10)))) %>%
    filter(age_cat == cat) %>%
    mutate(weight = full_sample_2_year_mec_exam_weight / 2,
           weight_norm = weight / mean(weight)) %>%
    ungroup()  %>%
    survey::svydesign(
      id = ~ masked_variance_pseudo_psu,
      strata = ~ masked_variance_pseudo_stratum,
      weights = ~ weight_norm,
      data = .,
      nest = TRUE
    )
  formula = as.formula(paste(var,  "~ gender"))
  test = svyttest(formula, survey_design)
  pvals = append(pvals, unname(test$p.value))
}

small %>% 
  select(metric, age_cat, Female = val_female, Male = val_male, mean_Female, mean_Male) %>% 
  bind_cols(p = pvals) %>% 
  mutate(pval = if_else(p < 0.001, "<0.001", format.pval(p, digits = 2))) %>% 
  # mutate(metric = factor(metric, 
  #                        labels = c("MVPA (min)", 
  #                                   "VPA (min)", "Time <1000",
  #                                   "Time [1000,2000]", "Time (2000,3000]", "Time (3000,4000]", "Time (4000,5000]", "Time (5000,6000]", "Time (6000,7000]", "Time (7000,8000]", "Time (8000,9000]", "Time >9000", "Total AC", "Mean AC", "Total MIMS"))) %>% 
  group_by(metric) %>% 
  gt() %>% 
  cols_hide(c("mean_Female", "mean_Male", "p")) %>% 
  tab_style(
    # Style to apply when conditions are met
    style = cell_fill(color = "lightpink"),
    locations = cells_body(
      rows = mean_Female > mean_Male & p < 0.05
    )
  ) %>% 
  tab_style(
    # Style to apply when conditions are met
    style = cell_fill(color = "lightblue"),
    locations = cells_body(
      rows = mean_Female < mean_Male & p < 0.05
    )
  )

age_cat	Female	Male	pval
total
(18,30]	2822818 (2775629,2870006)	2722939 (2675751,2809921)	0.01498
(30,40]	2911286 (2840893,2981680)	2690187 (2619794,2746393)	<0.001
(40,50]	2814679 (2736284,2893075)	2612294 (2533899,2686557)	<0.001
(50,60]	2714765 (2638819,2790712)	2437418 (2361472,2507407)	<0.001
(60,70]	2465049 (2401490,2528608)	2103828 (2040269,2166496)	<0.001
(70,80]	1905875 (1840890,1970860)	1706065 (1641080,1761126)	<0.001
mean
(18,30]	1963 (1930,1995)	1895 (1862,1955)	0.01498
(30,40]	2024 (1975,2074)	1871 (1822,1910)	<0.001
(40,50]	1957 (1902,2011)	1816 (1762,1868)	<0.001
(50,60]	1887 (1834,1940)	1695 (1642,1743)	<0.001
(60,70]	1713 (1669,1758)	1462 (1418,1506)	<0.001
(70,80]	1324 (1279,1370)	1185 (1140,1224)	<0.001
mvpa_montoye
(18,30]	332 (325,339)	315 (309,328)	0.01706
(30,40]	347 (338,356)	313 (303,321)	<0.001
(40,50]	336 (325,347)	305 (294,316)	<0.001
(50,60]	328 (316,340)	282 (270,292)	<0.001
(60,70]	295 (286,304)	238 (229,249)	<0.001
(70,80]	211 (201,221)	175 (165,184)	<0.001
mod_montoye
(18,30]	121 (118,124)	120 (117,124)	0.01706
(30,40]	129 (126,132)	124 (121,127)	<0.001
(40,50]	130 (127,134)	125 (122,129)	<0.001
(50,60]	134 (129,139)	124 (119,128)	<0.001
(60,70]	129 (125,132)	114 (110,118)	<0.001
(70,80]	104 (99,108)	91 (86,94)	<0.001
vig_montoye
(18,30]	173 (169,178)	158 (154,166)	0.00703
(30,40]	178 (171,185)	149 (142,156)	<0.001
(40,50]	166 (158,174)	139 (131,146)	<0.001
(50,60]	152 (143,160)	118 (109,125)	<0.001
(60,70]	125 (118,132)	85 (78,91)	<0.001
(70,80]	69 (64,74)	49 (44,54)	<0.001
mvpa_montoye_bout
(18,30]	122 (118,127)	114 (109,123)	0.00703
(30,40]	130 (122,137)	107 (100,113)	<0.001
(40,50]	121 (113,129)	101 (92,108)	<0.001
(50,60]	114 (104,123)	84 (74,91)	<0.001
(60,70]	93 (86,99)	60 (53,66)	<0.001
(70,80]	49 (45,53)	33 (28,37)	<0.001
vig_montoye_bout
(18,30]	45 (43,47)	40 (38,45)	0.64368
(30,40]	47 (43,51)	35 (31,38)	0.00996
(40,50]	41 (37,45)	32 (28,35)	0.05577
(50,60]	35 (30,40)	23 (18,27)	0.00356
(60,70]	25 (21,30)	13 (9,16)	<0.001
(70,80]	9 (7,10)	5 (4,7)	<0.001
x0_100
(18,30]	580 (573,588)	555 (547,565)	0.64368
(30,40]	557 (550,565)	541 (533,549)	0.00996
(40,50]	562 (552,573)	549 (538,558)	0.05577
(50,60]	564 (554,574)	558 (548,570)	0.00356
(60,70]	587 (578,597)	599 (589,609)	<0.001
(70,80]	644 (632,656)	647 (635,658)	<0.001
x100_500
(18,30]	113 (110,115)	129 (126,133)	<0.001
(30,40]	113 (111,116)	132 (130,135)	<0.001
(40,50]	116 (113,120)	132 (128,136)	<0.001
(50,60]	116 (113,120)	137 (133,140)	<0.001
(60,70]	127 (124,130)	147 (144,152)	<0.001
(70,80]	147 (143,152)	162 (158,165)	<0.001
x500_1e_03
(18,30]	107 (104,110)	120 (118,123)	<0.001
(30,40]	108 (106,110)	125 (123,128)	<0.001
(40,50]	110 (107,113)	125 (123,129)	<0.001
(50,60]	110 (107,113)	130 (126,132)	<0.001
(60,70]	113 (111,115)	133 (130,137)	<0.001
(70,80]	121 (118,123)	134 (132,137)	<0.001
x1e_03_2e_03
(18,30]	152 (149,155)	163 (160,166)	0.02708
(30,40]	154 (152,157)	167 (165,170)	<0.001
(40,50]	155 (152,158)	167 (164,171)	<0.001
(50,60]	157 (153,161)	169 (165,174)	<0.001
(60,70]	155 (151,158)	165 (161,170)	<0.001
(70,80]	158 (155,160)	164 (161,168)	<0.001
x2e_03_3e_03
(18,30]	109 (107,110)	110 (109,113)	0.02708
(30,40]	111 (110,113)	113 (112,116)	<0.001
(40,50]	111 (109,114)	114 (111,116)	<0.001
(50,60]	114 (112,117)	115 (113,118)	<0.001
(60,70]	113 (110,116)	111 (108,113)	<0.001
(70,80]	112 (110,115)	112 (109,115)	<0.001
x3e_03_4e_03
(18,30]	87 (86,89)	87 (85,89)	0.01478
(30,40]	92 (90,93)	91 (89,93)	<0.001
(40,50]	92 (90,94)	92 (90,94)	0.00171
(50,60]	96 (94,99)	93 (90,96)	<0.001
(60,70]	95 (93,98)	91 (88,93)	<0.001
(70,80]	90 (87,92)	85 (82,88)	<0.001
x4e_03_5e_03
(18,30]	75 (74,77)	75 (73,78)	0.01478
(30,40]	80 (79,82)	78 (76,80)	<0.001
(40,50]	81 (79,83)	79 (77,81)	0.00171
(50,60]	84 (81,87)	79 (76,81)	<0.001
(60,70]	82 (80,84)	74 (72,77)	<0.001
(70,80]	68 (65,71)	61 (58,63)	<0.001
x5e_03_6e_03
(18,30]	65 (63,66)	63 (62,66)	<0.001
(30,40]	69 (67,71)	65 (63,66)	0.00584
(40,50]	69 (67,71)	64 (62,67)	0.06039
(50,60]	70 (67,73)	62 (59,65)	0.34903
(60,70]	64 (63,66)	53 (51,56)	0.08821
(70,80]	46 (44,48)	37 (35,40)	0.64263
x6e_03_7e_03
(18,30]	52 (50,53)	49 (48,51)	<0.001
(30,40]	55 (53,56)	49 (47,51)	0.00584
(40,50]	53 (51,55)	47 (45,49)	0.06039
(50,60]	51 (48,54)	42 (39,44)	0.34903
(60,70]	44 (42,46)	33 (31,35)	0.08821
(70,80]	27 (25,29)	20 (18,22)	0.64263
x7e_03_8e_03
(18,30]	38 (37,39)	34 (33,36)	<0.001
(30,40]	39 (37,40)	32 (31,34)	<0.001
(40,50]	36 (35,38)	30 (28,32)	<0.001
(50,60]	33 (31,35)	25 (23,27)	<0.001
(60,70]	27 (25,29)	18 (16,19)	<0.001
(70,80]	14 (13,16)	9 (8,10)	<0.001
x8e_03_9e_03
(18,30]	25 (24,26)	22 (21,23)	<0.001
(30,40]	25 (24,26)	20 (18,21)	<0.001
(40,50]	22 (21,24)	17 (16,19)	<0.001
(50,60]	19 (18,21)	14 (12,15)	<0.001
(60,70]	15 (14,16)	8 (7,9)	<0.001
(70,80]	7 (6,8)	4 (3,5)	<0.001
x9e_03_inf
(18,30]	35 (34,37)	30 (28,32)	<0.001
(30,40]	35 (31,38)	25 (22,27)	<0.001
(40,50]	29 (27,32)	22 (20,25)	<0.001
(50,60]	23 (21,25)	16 (14,18)	<0.001
(60,70]	16 (14,18)	9 (6,10)	<0.001
(70,80]	5 (5,6)	4 (3,4)	<0.001

Stratified by wave

2011-2012

Code

survey_design =
  df_small_gh %>%
  filter(data_release_cycle == 7) %>%
  mutate(age_cat = cut(age_in_years_at_screening,
                       breaks=c(18, 30, seq(40, 80, 10)))) %>%
  mutate(weight = full_sample_2_year_mec_exam_weight,
         weight_norm = weight / mean(weight)) %>%
  ungroup()  %>%
  survey::svydesign(
    id = ~ masked_variance_pseudo_psu,
    strata = ~ masked_variance_pseudo_stratum,
    weights = ~ weight_norm,
    data = .,
    nest = TRUE
  )

calc_by_age =
  function(var, df) {
    formula = as.formula(paste("~", var))
    total_by_age_gender = svyby(formula,
                                ~ age_cat + gender,
                                df,
                                svymean) %>%
      rename(mean = contains(var)) %>%
      mutate(metric = var)
  }

means_df =
  map_dfr(.x = pa %>% select(total, 
         mean,
         mvpa_montoye, 
         mod_montoye,
         vig_montoye,
         mvpa_montoye_bout,
         vig_montoye_bout,
         starts_with("x")) %>% colnames(),
          .f = calc_by_age, df = survey_design)
small = 
  means_df %>%
  pivot_wider(names_from = gender, values_from = c(mean, se)) %>% 
  mutate(val_female = paste0(round(mean_Female), " (", round(mean_Female - 1.96 * se_Female), ",", round(mean_Female + 1.96 * se_Female), ")"),
         val_male = paste0(round(mean_Male), " (", round(mean_Male - 1.96 * se_Female), ",", round(mean_Male + 1.96 * se_Male), ")"))

# test = svyttest(mvpa_montoye ~ gender, survey_design)
pvals = c()
for(i in 1:nrow(small)){
  row = means_df[i,]
  cat = row$age_cat
  var = row$metric
  survey_design =
    df_small_gh %>%
    filter(age_in_years_at_screening >= 18 & data_release_cycle == 7) %>%
    mutate(age_cat = cut(age_in_years_at_screening,
                         breaks=c(18, 30, seq(40, 80, 10)))) %>%
    filter(age_cat == cat) %>%
    mutate(weight = full_sample_2_year_mec_exam_weight / 2,
           weight_norm = weight / mean(weight)) %>%
    ungroup()  %>%
    survey::svydesign(
      id = ~ masked_variance_pseudo_psu,
      strata = ~ masked_variance_pseudo_stratum,
      weights = ~ weight_norm,
      data = .,
      nest = TRUE
    )
  formula = as.formula(paste(var,  "~ gender"))
  test = svyttest(formula, survey_design)
  pvals = append(pvals, unname(test$p.value))
}

small %>% 
  select(metric, age_cat, Female = val_female, Male = val_male, mean_Female, mean_Male) %>% 
  bind_cols(p = pvals) %>% 
  mutate(pval = if_else(p < 0.001, "<0.001", format.pval(p, digits = 2))) %>% 
  # mutate(metric = factor(metric, 
  #                        labels = c("MVPA (min)", 
  #                                   "VPA (min)", "Time <1000",
  #                                   "Time [1000,2000]", "Time (2000,3000]", "Time (3000,4000]", "Time (4000,5000]", "Time (5000,6000]", "Time (6000,7000]", "Time (7000,8000]", "Time (8000,9000]", "Time >9000", "Total AC", "Mean AC", "Total MIMS"))) %>% 
  group_by(metric) %>% 
  gt() %>% 
  cols_hide(c("mean_Female", "mean_Male", "p")) %>% 
  tab_style(
    # Style to apply when conditions are met
    style = cell_fill(color = "lightpink"),
    locations = cells_body(
      rows = mean_Female > mean_Male & p < 0.05
    )
  ) %>% 
  tab_style(
    # Style to apply when conditions are met
    style = cell_fill(color = "lightblue"),
    locations = cells_body(
      rows = mean_Female < mean_Male & p < 0.05
    )
  )

age_cat	Female	Male	pval
total
(18,30]	2760827 (2703034,2818621)	2704620 (2646826,2838398)	0.38584
(30,40]	2928179 (2838858,3017500)	2738119 (2648798,2817110)	0.00115
(40,50]	2842576 (2713747,2971406)	2616878 (2488049,2736548)	0.03035
(50,60]	2794832 (2673840,2915825)	2419536 (2298543,2527517)	<0.001
(60,70]	2440773 (2360832,2520714)	2144570 (2064629,2209794)	<0.001
(70,80]	1890705 (1820594,1960816)	1675078 (1604967,1755619)	<0.001
mean
(18,30]	1920 (1880,1960)	1882 (1842,1975)	0.38584
(30,40]	2036 (1973,2099)	1904 (1842,1960)	0.00115
(40,50]	1976 (1886,2065)	1819 (1730,1902)	0.03035
(50,60]	1943 (1859,2027)	1682 (1598,1757)	<0.001
(60,70]	1696 (1640,1751)	1490 (1435,1536)	<0.001
(70,80]	1314 (1265,1362)	1164 (1115,1220)	<0.001
mvpa_montoye
(18,30]	323 (314,332)	310 (301,330)	0.40087
(30,40]	350 (338,363)	318 (305,329)	0.00120
(40,50]	341 (322,359)	306 (288,324)	0.03035
(50,60]	338 (319,357)	279 (260,295)	<0.001
(60,70]	292 (281,304)	244 (232,253)	<0.001
(70,80]	210 (199,220)	169 (158,181)	<0.001
mod_montoye
(18,30]	120 (116,125)	118 (113,125)	0.40087
(30,40]	131 (127,134)	125 (121,128)	0.00120
(40,50]	132 (127,138)	126 (121,133)	0.03035
(50,60]	134 (128,141)	124 (118,130)	<0.001
(60,70]	130 (126,135)	116 (111,120)	<0.001
(70,80]	104 (98,110)	88 (81,93)	<0.001
vig_montoye
(18,30]	165 (159,171)	155 (149,167)	0.20492
(30,40]	179 (170,188)	154 (145,162)	<0.001
(40,50]	167 (155,180)	140 (127,151)	0.01984
(50,60]	162 (148,176)	114 (100,125)	<0.001
(60,70]	120 (111,129)	88 (79,92)	<0.001
(70,80]	67 (63,71)	46 (42,53)	<0.001
mvpa_montoye_bout
(18,30]	114 (106,121)	109 (102,123)	0.20492
(30,40]	132 (123,141)	111 (101,118)	<0.001
(40,50]	121 (108,134)	101 (88,112)	0.01984
(50,60]	126 (110,142)	80 (65,91)	<0.001
(60,70]	87 (79,96)	61 (53,67)	<0.001
(70,80]	48 (44,52)	31 (27,37)	<0.001
vig_montoye_bout
(18,30]	41 (38,44)	39 (36,45)	0.59327
(30,40]	48 (43,52)	38 (33,42)	0.03672
(40,50]	40 (34,47)	31 (25,36)	0.16619
(50,60]	42 (33,50)	22 (13,26)	0.02537
(60,70]	22 (17,27)	13 (9,15)	<0.001
(70,80]	8 (7,9)	5 (5,8)	<0.001
x0_100
(18,30]	582 (570,593)	553 (541,570)	0.59327
(30,40]	552 (544,561)	536 (527,549)	0.03672
(40,50]	553 (539,567)	549 (535,565)	0.16619
(50,60]	559 (544,574)	557 (542,573)	0.02537
(60,70]	582 (568,596)	586 (572,603)	<0.001
(70,80]	646 (630,662)	647 (631,663)	<0.001
x100_500
(18,30]	114 (109,118)	130 (126,137)	0.08441
(30,40]	114 (110,117)	129 (126,133)	<0.001
(40,50]	116 (110,122)	129 (123,136)	0.00704
(50,60]	114 (110,119)	138 (133,144)	<0.001
(60,70]	129 (124,134)	149 (145,153)	<0.001
(70,80]	149 (142,156)	162 (155,168)	<0.001
x500_1e_03
(18,30]	109 (104,114)	122 (117,126)	0.08441
(30,40]	109 (106,111)	125 (122,128)	<0.001
(40,50]	110 (106,114)	124 (120,128)	0.00704
(50,60]	108 (104,112)	130 (126,135)	<0.001
(60,70]	115 (112,119)	135 (132,140)	<0.001
(70,80]	120 (116,124)	137 (134,141)	<0.001
x1e_03_2e_03
(18,30]	155 (150,160)	165 (161,169)	0.47957
(30,40]	154 (151,158)	169 (166,173)	0.00310
(40,50]	157 (154,161)	168 (164,173)	0.02780
(50,60]	156 (150,163)	171 (165,180)	<0.001
(60,70]	156 (150,163)	167 (161,175)	<0.001
(70,80]	156 (151,160)	168 (164,173)	<0.001
x2e_03_3e_03
(18,30]	110 (108,112)	111 (109,115)	0.47957
(30,40]	112 (110,114)	114 (112,117)	0.00310
(40,50]	114 (110,117)	114 (111,119)	0.02780
(50,60]	114 (110,117)	116 (113,121)	<0.001
(60,70]	115 (111,119)	111 (106,115)	<0.001
(70,80]	112 (108,116)	112 (108,117)	<0.001
x3e_03_4e_03
(18,30]	88 (86,90)	87 (84,91)	0.54589
(30,40]	92 (90,95)	91 (89,93)	<0.001
(40,50]	94 (91,98)	93 (89,97)	0.03529
(50,60]	95 (92,99)	94 (90,98)	0.00107
(60,70]	97 (94,101)	91 (88,95)	0.00433
(70,80]	90 (86,94)	84 (79,88)	0.03401
x4e_03_5e_03
(18,30]	75 (72,78)	74 (71,78)	0.54589
(30,40]	81 (79,84)	78 (76,80)	<0.001
(40,50]	83 (79,86)	79 (76,83)	0.03529
(50,60]	84 (80,88)	79 (75,83)	0.00107
(60,70]	83 (80,87)	75 (72,78)	0.00433
(70,80]	68 (64,73)	59 (55,63)	0.03401
x5e_03_6e_03
(18,30]	64 (61,66)	62 (59,66)	0.00666
(30,40]	70 (68,72)	65 (63,67)	0.01437
(40,50]	70 (66,73)	65 (61,69)	0.73753
(50,60]	71 (66,75)	61 (57,65)	0.84044
(60,70]	64 (62,66)	55 (52,57)	0.73466
(70,80]	46 (43,49)	35 (33,39)	0.89408
x6e_03_7e_03
(18,30]	50 (48,52)	48 (46,51)	0.00666
(30,40]	55 (53,58)	50 (47,52)	0.01437
(40,50]	53 (50,57)	47 (44,51)	0.73753
(50,60]	54 (49,58)	41 (37,44)	0.84044
(60,70]	43 (41,46)	33 (31,35)	0.73466
(70,80]	27 (25,28)	18 (17,21)	0.89408
x7e_03_8e_03
(18,30]	36 (35,38)	33 (32,36)	0.00119
(30,40]	39 (37,41)	33 (31,35)	<0.001
(40,50]	37 (34,40)	30 (27,33)	0.00828
(50,60]	36 (32,39)	24 (21,27)	<0.001
(60,70]	26 (24,28)	18 (15,19)	<0.001
(70,80]	14 (13,15)	9 (8,10)	0.00982
x8e_03_9e_03
(18,30]	23 (22,25)	21 (20,23)	0.00119
(30,40]	25 (24,26)	20 (19,22)	<0.001
(40,50]	23 (20,25)	18 (15,19)	0.00828
(50,60]	21 (19,24)	13 (11,15)	<0.001
(60,70]	14 (12,16)	8 (7,9)	<0.001
(70,80]	7 (6,7)	4 (3,5)	0.00982
x9e_03_inf
(18,30]	33 (30,35)	30 (28,33)	<0.001
(30,40]	34 (30,38)	27 (24,30)	<0.001
(40,50]	29 (26,32)	22 (18,25)	<0.001
(50,60]	26 (22,29)	15 (11,18)	<0.001
(60,70]	14 (11,17)	10 (7,11)	<0.001
(70,80]	5 (4,5)	4 (3,5)	<0.001

2013-2014

Code

survey_design =
  df_small_gh %>%
  filter(data_release_cycle == 8) %>%
  mutate(age_cat = cut(age_in_years_at_screening,
                       breaks=c(18, 30, seq(40, 80, 10)))) %>%
  mutate(weight = full_sample_2_year_mec_exam_weight,
         weight_norm = weight / mean(weight)) %>%
  ungroup()  %>%
  survey::svydesign(
    id = ~ masked_variance_pseudo_psu,
    strata = ~ masked_variance_pseudo_stratum,
    weights = ~ weight_norm,
    data = .,
    nest = TRUE
  )

calc_by_age =
  function(var, df) {
    formula = as.formula(paste("~", var))
    total_by_age_gender = svyby(formula,
                                ~ age_cat + gender,
                                df,
                                svymean) %>%
      rename(mean = contains(var)) %>%
      mutate(metric = var)
  }

means_df =
  map_dfr(.x = pa %>% select(total, 
         mean,
         mvpa_montoye, 
         mod_montoye,
         vig_montoye,
         mvpa_montoye_bout,
         vig_montoye_bout,
         starts_with("x")) %>% colnames(),
          .f = calc_by_age, df = survey_design)

small = 
  means_df %>%
  pivot_wider(names_from = gender, values_from = c(mean, se)) %>% 
  mutate(val_female = paste0(round(mean_Female), " (", round(mean_Female - 1.96 * se_Female), ",", round(mean_Female + 1.96 * se_Female), ")"),
         val_male = paste0(round(mean_Male), " (", round(mean_Male - 1.96 * se_Female), ",", round(mean_Male + 1.96 * se_Male), ")"))

# test = svyttest(mvpa_montoye ~ gender, survey_design)
pvals = c()
for(i in 1:nrow(small)){
  row = means_df[i,]
  cat = row$age_cat
  var = row$metric
  survey_design =
    df_small_gh %>%
    filter(age_in_years_at_screening >= 18 & data_release_cycle == 8) %>%
    mutate(age_cat = cut(age_in_years_at_screening,
                         breaks=c(18, 30, seq(40, 80, 10)))) %>%
    filter(age_cat == cat) %>%
    mutate(weight = full_sample_2_year_mec_exam_weight / 2,
           weight_norm = weight / mean(weight)) %>%
    ungroup()  %>%
    survey::svydesign(
      id = ~ masked_variance_pseudo_psu,
      strata = ~ masked_variance_pseudo_stratum,
      weights = ~ weight_norm,
      data = .,
      nest = TRUE
    )
  formula = as.formula(paste(var,  "~ gender"))
  test = svyttest(formula, survey_design)
  pvals = append(pvals, unname(test$p.value))
}

small %>% 
  select(metric, age_cat, Female = val_female, Male = val_male, mean_Female, mean_Male) %>% 
  bind_cols(p = pvals) %>% 
  mutate(pval = if_else(p < 0.001, "<0.001", format.pval(p, digits = 2))) %>% 
  # mutate(metric = factor(metric, 
  #                        labels = c("MVPA (min)", 
  #                                   "VPA (min)", "Time <1000",
  #                                   "Time [1000,2000]", "Time (2000,3000]", "Time (3000,4000]", "Time (4000,5000]", "Time (5000,6000]", "Time (6000,7000]", "Time (7000,8000]", "Time (8000,9000]", "Time >9000", "Total AC", "Mean AC", "Total MIMS"))) %>% 
  group_by(metric) %>% 
  gt() %>% 
  cols_hide(c("mean_Female", "mean_Male", "p")) %>% 
  tab_style(
    # Style to apply when conditions are met
    style = cell_fill(color = "lightpink"),
    locations = cells_body(
      rows = mean_Female > mean_Male & p < 0.05
    )
  ) %>% 
  tab_style(
    # Style to apply when conditions are met
    style = cell_fill(color = "lightblue"),
    locations = cells_body(
      rows = mean_Female < mean_Male & p < 0.05
    )
  )

age_cat	Female	Male	pval
total
(18,30]	2882509 (2808573,2956446)	2742097 (2668160,2853457)	0.00689
(30,40]	2895061 (2787659,3002462)	2639481 (2532080,2719143)	0.00645
(40,50]	2786316 (2696997,2875635)	2607407 (2518088,2692178)	0.00323
(50,60]	2631658 (2553190,2710126)	2457263 (2378795,2543668)	0.00888
(60,70]	2489119 (2393179,2585059)	2065497 (1969557,2170579)	<0.001
(70,80]	1920711 (1811994,2029429)	1734750 (1626032,1810594)	0.00703
mean
(18,30]	2004 (1953,2056)	1909 (1857,1986)	0.00689
(30,40]	2013 (1939,2088)	1835 (1761,1891)	0.00645
(40,50]	1938 (1876,2000)	1813 (1751,1872)	0.00323
(50,60]	1829 (1774,1883)	1709 (1654,1768)	0.00888
(60,70]	1730 (1663,1797)	1435 (1368,1508)	<0.001
(70,80]	1335 (1259,1411)	1206 (1130,1258)	0.00703
mvpa_montoye
(18,30]	341 (331,350)	321 (311,336)	0.00761
(30,40]	344 (330,358)	308 (294,319)	0.00630
(40,50]	332 (319,344)	303 (291,316)	0.00311
(50,60]	317 (303,331)	286 (272,299)	0.00916
(60,70]	298 (285,311)	233 (220,252)	<0.001
(70,80]	212 (195,229)	182 (165,194)	0.00714
mod_montoye
(18,30]	122 (119,125)	122 (119,127)	0.00761
(30,40]	128 (124,132)	124 (120,128)	0.00630
(40,50]	128 (124,132)	124 (121,129)	0.00311
(50,60]	133 (127,140)	124 (118,130)	0.00916
(60,70]	127 (123,131)	111 (107,119)	<0.001
(70,80]	103 (97,110)	93 (87,99)	0.00714
vig_montoye
(18,30]	181 (175,188)	162 (155,173)	0.00555
(30,40]	176 (165,188)	145 (133,154)	0.00343
(40,50]	165 (156,174)	139 (130,148)	<0.001
(50,60]	141 (134,149)	122 (115,130)	0.00541
(60,70]	130 (120,141)	83 (73,94)	<0.001
(70,80]	71 (61,80)	52 (42,58)	0.00887
mvpa_montoye_bout
(18,30]	131 (125,137)	118 (113,130)	0.00555
(30,40]	127 (117,138)	102 (92,112)	0.00343
(40,50]	121 (111,131)	100 (91,111)	<0.001
(50,60]	101 (93,109)	88 (79,96)	0.00541
(60,70]	98 (88,108)	59 (48,69)	<0.001
(70,80]	50 (42,57)	34 (27,40)	0.00887
vig_montoye_bout
(18,30]	50 (46,53)	42 (38,48)	0.93807
(30,40]	47 (40,53)	32 (26,37)	0.11891
(40,50]	41 (36,47)	33 (27,38)	0.15608
(50,60]	28 (25,31)	25 (22,31)	0.06455
(60,70]	29 (23,35)	13 (7,17)	0.00153
(70,80]	10 (7,12)	5 (3,7)	0.03433
x0_100
(18,30]	579 (569,589)	556 (547,568)	0.93807
(30,40]	562 (550,574)	546 (534,556)	0.11891
(40,50]	572 (556,588)	549 (533,558)	0.15608
(50,60]	570 (557,582)	559 (546,576)	0.06455
(60,70]	593 (580,607)	611 (598,623)	0.00153
(70,80]	642 (625,659)	647 (630,662)	0.03433
x100_500
(18,30]	112 (110,115)	127 (124,132)	<0.001
(30,40]	112 (109,116)	135 (131,139)	0.00217
(40,50]	117 (113,120)	134 (131,139)	<0.001
(50,60]	119 (113,124)	136 (131,140)	0.00410
(60,70]	125 (121,129)	145 (141,153)	<0.001
(70,80]	146 (141,150)	161 (157,164)	0.00346
x500_1e_03
(18,30]	105 (104,106)	119 (118,122)	<0.001
(30,40]	108 (105,111)	126 (123,129)	0.00217
(40,50]	110 (107,113)	127 (123,131)	<0.001
(50,60]	111 (107,116)	129 (124,132)	0.00410
(60,70]	111 (108,113)	130 (128,136)	<0.001
(70,80]	121 (118,124)	132 (128,137)	0.00346
x1e_03_2e_03
(18,30]	149 (148,151)	160 (159,164)	0.01716
(30,40]	154 (151,158)	165 (162,168)	0.00848
(40,50]	153 (148,158)	166 (162,171)	0.00181
(50,60]	158 (153,163)	167 (162,170)	0.02001
(60,70]	153 (149,157)	162 (158,168)	<0.001
(70,80]	159 (156,163)	161 (157,165)	0.00329
x2e_03_3e_03
(18,30]	107 (105,109)	109 (107,112)	0.01716
(30,40]	111 (108,114)	113 (110,116)	0.00848
(40,50]	109 (106,112)	113 (110,115)	0.00181
(50,60]	115 (112,118)	114 (111,117)	0.02001
(60,70]	111 (107,115)	111 (107,114)	<0.001
(70,80]	113 (109,116)	111 (108,115)	0.00329
x3e_03_4e_03
(18,30]	87 (85,88)	87 (85,90)	0.00961
(30,40]	91 (88,94)	91 (88,93)	0.00471
(40,50]	90 (88,92)	91 (89,94)	0.02262
(50,60]	97 (93,101)	92 (88,96)	0.42677
(60,70]	94 (90,97)	90 (87,93)	<0.001
(70,80]	89 (85,93)	87 (82,90)	0.00839
x4e_03_5e_03
(18,30]	76 (74,77)	76 (74,79)	0.00961
(30,40]	80 (77,82)	78 (75,80)	0.00471
(40,50]	80 (78,82)	78 (76,81)	0.02262
(50,60]	84 (80,88)	78 (74,82)	0.42677
(60,70]	80 (77,83)	73 (70,78)	<0.001
(70,80]	68 (64,72)	63 (59,66)	0.00839
x5e_03_6e_03
(18,30]	66 (64,68)	65 (63,67)	0.00350
(30,40]	68 (66,70)	64 (62,66)	0.10196
(40,50]	68 (65,70)	64 (61,66)	0.01199
(50,60]	68 (64,72)	63 (59,66)	0.22575
(60,70]	64 (62,67)	52 (49,56)	0.01183
(70,80]	46 (42,50)	39 (35,42)	0.60374
x6e_03_7e_03
(18,30]	54 (52,55)	51 (49,54)	0.00350
(30,40]	54 (51,56)	48 (46,51)	0.10196
(40,50]	52 (50,55)	46 (44,49)	0.01199
(50,60]	49 (46,52)	43 (40,45)	0.22575
(60,70]	45 (42,47)	32 (29,36)	0.01183
(70,80]	27 (24,30)	21 (18,23)	0.60374
x7e_03_8e_03
(18,30]	40 (38,41)	36 (34,38)	<0.001
(30,40]	38 (36,40)	32 (30,34)	<0.001
(40,50]	36 (34,38)	29 (27,32)	<0.001
(50,60]	31 (29,33)	26 (24,28)	<0.001
(60,70]	28 (26,31)	17 (15,20)	<0.001
(70,80]	15 (13,17)	10 (8,11)	<0.001
x8e_03_9e_03
(18,30]	26 (25,28)	22 (21,24)	<0.001
(30,40]	25 (23,27)	19 (17,21)	<0.001
(40,50]	22 (21,24)	17 (16,19)	<0.001
(50,60]	18 (16,19)	14 (13,16)	<0.001
(60,70]	16 (14,18)	8 (6,10)	<0.001
(70,80]	7 (6,9)	4 (3,5)	<0.001
x9e_03_inf
(18,30]	38 (35,40)	30 (27,33)	<0.001
(30,40]	35 (30,40)	23 (18,26)	<0.001
(40,50]	30 (27,33)	23 (20,27)	<0.001
(50,60]	20 (18,21)	16 (15,19)	<0.001
(60,70]	18 (14,22)	8 (4,9)	<0.001
(70,80]	6 (5,8)	4 (2,4)	0.00206

Plots

Distributions

Code

pa %>% 
  left_join(demo_gh) %>% 
  filter(age_in_years_at_screening >= 18) %>% 
  pivot_longer(cols = c(contains("montoye"), mean)) %>% 
  mutate(name = factor(name, labels = c("LPA", "Mean AC", "MPA", "MVPA", "MVPA 10 min bout", "Sedentary", "VPA", "VPA 10 min bout"))) %>% 
  ggplot(aes(x = value, color = gender)) + 
  geom_density() + 
  facet_wrap(.~name, scales = "free") + 
  scale_color_manual(values = c(col1, col2), name = "")+
  theme(legend.position = c(0.9, 0.2))

Smoothed, unweighted

Code

pa %>%
  left_join(demo_gh) %>% 
  filter(age_in_years_at_screening >= 18) %>% 
  pivot_longer(cols = contains("bout")) %>% 
  mutate(name = factor(name, labels = c("MVPA 10 min Bout", "VPA 10 min Bout"))) %>% 
  ggplot(aes(x = age_in_years_at_screening, y = value, color = gender)) +
  geom_jitter(alpha = .1, size = .7, width = .25) +
  geom_smooth(se = F) +
  facet_grid(name~., scales = "free_y") +
  labs(x = "Age", y = "Minutes", color = "")+
  scale_color_manual(values = c(col1, col2)) +
  theme(legend.position = c(0.9, 0.9),
        legend.title = element_blank())

Code

pa %>%
  left_join(demo_gh) %>% 
  filter(age_in_years_at_screening >= 18) %>% 
  pivot_longer(cols = contains("montoye") & !contains("bout")) %>% 
  mutate(name = factor(name, labels = c("Light", "Moderate", "MVPA", "Sedentary", "Vigorous"))) %>%
  ggplot(aes(x = age_in_years_at_screening, y = value, color = gender)) +
  geom_jitter(alpha = .1, size = .7, width = .25) +
  geom_smooth(se = F) +
  facet_grid(name~., scales = "free_y") + 
  labs(x = "Age", y = "Minutes", color = "")+
  scale_color_manual(values = c(col1, col2)) +
  theme(legend.position = c(0.9, 0.9),
        legend.title = element_blank())

Code

pa %>%
  left_join(demo_gh) %>% 
  filter(age_in_years_at_screening >= 18) %>% 
  pivot_longer(cols = starts_with("x")) %>% 
  ggplot(aes(x = age_in_years_at_screening, y = value, color = gender)) +
  geom_jitter(alpha = .1, size = .7, width = .25) +
  geom_smooth(se = F) +
  facet_wrap(name~., scales = "free_y") + 
  labs(x = "Age", y = "Minutes", color = "")+
  scale_color_manual(values = c(col1, col2)) +
  theme(legend.position = "bottom",
        legend.title = element_blank())

Smoothed, weighted

Code

survey_design =
  df_small_gh %>%
  filter(age_in_years_at_screening >= 18) %>%
  mutate(weight = full_sample_2_year_mec_exam_weight / 2,
         weight_norm = weight / mean(weight)) %>%
  ungroup()  %>%
  survey::svydesign(
    id = ~ masked_variance_pseudo_psu,
    strata = ~ masked_variance_pseudo_stratum,
    weights = ~ weight_norm,
    data = .,
    nest = TRUE
  )


calc_by_age =
  function(var, df) {
    # var = "total_oaksteps"
    formula = as.formula(paste("~", var))
    total_by_age_gender = svyby(formula,
                                ~ age_in_years_at_screening + gender,
                                df,
                                svymean) %>%
      rename(mean = contains(var)) %>%
      mutate(metric = var)
  }

means_df =
  map_dfr(.x = pa %>% select(total, 
         mean,
         mvpa_montoye, 
         mod_montoye,
         vig_montoye,
         mvpa_montoye_bout,
         vig_montoye_bout,
         starts_with("x")) %>% colnames(),
          .f = calc_by_age, df = survey_design)


# get smooths for means and confidence intervals
models = means_df %>%
  mutate(lb = mean - 1.96 * se,
         ub = mean + 1.96 * se) %>%
  tidyr::nest(data = -c(metric, gender)) %>%
  dplyr::mutate(
    # Perform loess calculation on each group
    m = purrr::map(data, loess,
                   formula = mean ~ age_in_years_at_screening, span = .75),
    # Retrieve the fitted values from each model
    fitted_mean = purrr::map(m, `[[`, "fitted"),
    l = purrr::map(data, loess,
                   formula = lb ~ age_in_years_at_screening, span = .75),
    # Retrieve the fitted values from each model
    fitted_lb = purrr::map(l, `[[`, "fitted"),
    u = purrr::map(data, loess,
                   formula = ub ~ age_in_years_at_screening, span = .75),
    # Retrieve the fitted values from each model
    fitted_ub = purrr::map(u, `[[`, "fitted")
  )

# Apply fitted y's as a new column
results = models %>%
  dplyr::select(-m, -l, -u) %>%
  tidyr::unnest(cols = c(data, contains("fitted")))

results %>%
  filter(!grepl("x", metric)) %>% 
  # mutate(metric = factor(metric, labels = c("MVPA", "VPA"))) %>% 
    ggplot(aes(x = age_in_years_at_screening, y = fitted_mean,
             ymin = fitted_lb, ymax = fitted_ub, color = gender, fill = gender)) +
  facet_wrap(. ~ metric, scales = "free_y") +
  geom_line(linewidth = 1) +
  geom_ribbon(alpha = .2, linetype = 0) +
  theme_light() +
  theme(legend.position = "bottom",
        legend.title = element_blank(),
        text =element_text(size = 14),
        strip.text = element_text(size = 12),
        panel.grid.minor = element_blank()) +
  labs(x = "Age (years)", y = "Mean Daily Minutes",
       title = "Smoothed survey weighted mean daily minutes by age ")+
  scale_x_continuous(breaks=seq(20,80,10)) +
  scale_y_continuous(breaks=seq(0,600,30)) + 
  scale_color_manual(values = c(col1, col2))+
  scale_fill_manual(values = c(col1, col2))

Code

results %>% 
  filter(grepl("x", metric)) %>% 
    ggplot(aes(x = age_in_years_at_screening, y = fitted_mean,
             ymin = fitted_lb, ymax = fitted_ub, color = gender, fill = gender)) +
  facet_wrap(. ~ metric, scales = "free_y") +
  geom_line(linewidth = 1) +
  geom_ribbon(alpha = .2, linetype = 0) +
  theme_light() +
  theme(legend.position = "bottom",
        legend.title = element_blank(),
        text =element_text(size = 14),
        strip.text = element_text(size = 12),
        panel.grid.minor = element_blank()) +
  labs(x = "Age (years)", y = "Mean Daily Minutes",
       title = "Smoothed survey weighted mean daily minutes by age")+
  scale_x_continuous(breaks=seq(20,80,10)) +
  scale_color_manual(values = c(col1, col2))+
  scale_fill_manual(values = c(col1, col2))

MIMS

Code

pa_mims = readRDS(here::here("data", "accelerometry", 
                              "summarized", "pa_summary_PAXMTSM.rds")) %>% 
  mutate(SEQN = as.character(SEQN)) %>% 
  group_by(SEQN) %>%
  summarize(across(-PAXDAYM, ~mean(.x, na.rm = TRUE))) 

Code

df_small_gh_mims = pa_mims %>% 
  left_join(demo_gh) %>% 
  filter(age_in_years_at_screening >= 18) %>%
  select(SEQN,
         data_release_cycle,
         gender,
         masked_variance_pseudo_psu,
         masked_variance_pseudo_stratum,
         full_sample_2_year_interview_weight,
         full_sample_2_year_mec_exam_weight,
         age_in_years_at_screening,
         total, 
         mean,
         starts_with("x"),
         missing_mims = na) 
labs = c("SEQN", "Data release cycle", "Gender",  
         "Pseudo PSU",
         "Psueudo stratum",
         "2 yr int weight", "2 yr exam weight",
         "Age", "Total MIMS", "Mean MIMS",
         "Time <2", "Time [2,5]","Time (5,10]","Time (10,15]", "Time (15,20]", "Time (20,30]", "Time (30,40]", "Time (40,50]", "Time (50,60]",  "Time >60", "Missing MIMS") 

names(labs) = colnames(df_small_gh_mims)

df_small_gh_mims =
  df_small_gh_mims %>%
  labelled::set_variable_labels(!!!labs)

df_svy =
  df_small_gh_mims %>% 
  mutate(WTMEC4YR = full_sample_2_year_mec_exam_weight/2,
         WTMEC4YR_norm = WTMEC4YR/mean(WTMEC4YR, na.rm = TRUE)) %>%
  svydesign(ids = ~masked_variance_pseudo_psu, weights = ~WTMEC4YR_norm,
            strata = ~masked_variance_pseudo_stratum, data=., nest=TRUE)

# survey weighted table
df_svy %>%
  tbl_svysummary(
    by = gender,
    include = -c(SEQN, data_release_cycle, masked_variance_pseudo_psu, masked_variance_pseudo_stratum, full_sample_2_year_interview_weight, full_sample_2_year_mec_exam_weight, WTMEC4YR, WTMEC4YR_norm),
    statistic = list(
      all_continuous() ~ "{mean} ({sd})",
      all_categorical() ~ "{n} ({p}%)"
    ),
    digits = all_continuous() ~ 0,
    missing_text = "Missing",
  ) %>%
  add_overall() %>%
  add_p()

Characteristic	Overall N = 9,298 1	Female N = 4,876 1	Male N = 4,422 1	p-value 2
Age	47 (17)	48 (18)	47 (17)	<0.001
Total MIMS	13,496 (3,819)	13,820 (3,766)	13,139 (3,845)	<0.001
Mean MIMS	9 (3)	10 (3)	9 (3)	<0.001
Time <2	592 (110)	596 (107)	587 (112)	0.003
Time [2,5]	151 (43)	142 (40)	160 (46)	<0.001
Time (5,10]	185 (42)	177 (39)	195 (44)	<0.001
Time (10,15]	134 (31)	132 (29)	136 (33)	0.002
Time (15,20]	110 (29)	110 (28)	109 (30)	0.065
Time (20,30]	167 (60)	173 (59)	161 (61)	<0.001
Time (30,40]	74 (47)	80 (48)	67 (45)	<0.001
Time (40,50]	19 (20)	21 (21)	17 (19)	<0.001
Time (50,60]	4 (7)	4 (7)	4 (7)	<0.001
Time >60	2 (6)	2 (5)	3 (8)	<0.001
Missing MIMS	20 (21)	22 (22)	18 (20)	0.10
Missing	6,805	3,734	3,072
1 Mean (SD)
2 Design-based KruskalWallis test

Code

survey_design =
  df_small_gh_mims %>%
  mutate(age_cat = cut(age_in_years_at_screening,
                       breaks=c(18, 30, seq(40, 80, 10)))) %>%
  mutate(weight = full_sample_2_year_mec_exam_weight / 2,
         weight_norm = weight / mean(weight)) %>%
  ungroup()  %>%
  survey::svydesign(
    id = ~ masked_variance_pseudo_psu,
    strata = ~ masked_variance_pseudo_stratum,
    weights = ~ weight_norm,
    data = .,
    nest = TRUE
  )

calc_by_age =
  function(var, df) {
    formula = as.formula(paste("~", var))
    total_by_age_gender = svyby(formula,
                                ~ age_cat + gender,
                                df,
                                svymean) %>%
      rename(mean = contains(var)) %>%
      mutate(metric = var)
  }

means_df =
  map_dfr(.x = pa_mims %>% select(starts_with("x"), total, mean) %>% colnames(),
          .f = calc_by_age, df = survey_design)


small = 
  means_df %>%
  pivot_wider(names_from = gender, values_from = c(mean, se)) %>% 
  mutate(val_female = paste0(round(mean_Female), " (", round(mean_Female - 1.96 * se_Female), ",", round(mean_Female + 1.96 * se_Female), ")"),
         val_male = paste0(round(mean_Male), " (", round(mean_Male - 1.96 * se_Female), ",", round(mean_Male + 1.96 * se_Male), ")"))

# test = svyttest(mvpa_montoye ~ gender, survey_design)
pvals = c()
for(i in 1:nrow(small)){
  row = means_df[i,]
  cat = row$age_cat
  var = row$metric
  survey_design =
    df_small_gh_mims %>%
    mutate(age_cat = cut(age_in_years_at_screening,
                         breaks=c(18, 30, seq(40, 80, 10)))) %>%
    filter(age_cat == cat) %>%
    mutate(weight = full_sample_2_year_mec_exam_weight / 2,
           weight_norm = weight / mean(weight)) %>%
    ungroup()  %>%
    survey::svydesign(
      id = ~ masked_variance_pseudo_psu,
      strata = ~ masked_variance_pseudo_stratum,
      weights = ~ weight_norm,
      data = .,
      nest = TRUE
    )
  formula = as.formula(paste(var,  "~ gender"))
  test = svyttest(formula, survey_design)
  pvals = append(pvals, unname(test$p.value))
}

small %>% 
  select(metric, age_cat, Female = val_female, Male = val_male, mean_Female, mean_Male) %>% 
  bind_cols(p = pvals) %>% 
  mutate(pval = if_else(p < 0.001, "<0.001", format.pval(p, digits = 2))) %>% 
  # mutate(metric = factor(metric, 
  #                        labels = c("MVPA (min)", 
  #                                   "VPA (min)", "Time <1000",
  #                                   "Time [1000,2000]", "Time (2000,3000]", "Time (3000,4000]", "Time (4000,5000]", "Time (5000,6000]", "Time (6000,7000]", "Time (7000,8000]", "Time (8000,9000]", "Time >9000", "Total AC", "Mean AC", "Total MIMS"))) %>% 
  group_by(metric) %>% 
  gt() %>% 
  cols_hide(c("mean_Female", "mean_Male", "p")) %>% 
  tab_style(
    # Style to apply when conditions are met
    style = cell_fill(color = "lightpink"),
    locations = cells_body(
      rows = mean_Female > mean_Male & p < 0.05
    )
  ) %>% 
  tab_style(
    # Style to apply when conditions are met
    style = cell_fill(color = "lightblue"),
    locations = cells_body(
      rows = mean_Female < mean_Male & p < 0.05
    )
  )

age_cat	Female	Male	pval
x0_2
(18,30]	594 (586,601)	569 (561,579)	<0.001
(30,40]	571 (563,579)	555 (547,563)	0.01062
(40,50]	578 (567,588)	567 (556,577)	0.13192
(50,60]	581 (571,591)	579 (569,592)	0.77927
(60,70]	607 (597,618)	623 (613,635)	0.02304
(70,80]	667 (655,678)	678 (666,689)	0.09110
x2_5
(18,30]	134 (131,137)	152 (149,157)	<0.001
(30,40]	135 (132,138)	156 (153,160)	0.01062
(40,50]	138 (135,142)	154 (151,160)	0.13192
(50,60]	138 (134,142)	161 (157,165)	0.77927
(60,70]	148 (145,152)	170 (166,175)	0.02304
(70,80]	167 (163,172)	181 (177,185)	0.09110
x5_10
(18,30]	173 (169,176)	190 (186,193)	<0.001
(30,40]	174 (172,177)	196 (194,200)	<0.001
(40,50]	177 (173,180)	194 (191,199)	<0.001
(50,60]	178 (173,183)	199 (194,205)	<0.001
(60,70]	177 (173,181)	198 (194,204)	<0.001
(70,80]	183 (179,186)	192 (189,196)	<0.001
x10_15
(18,30]	129 (127,131)	133 (131,136)	<0.001
(30,40]	133 (130,135)	138 (136,141)	<0.001
(40,50]	132 (130,135)	139 (136,143)	<0.001
(50,60]	135 (132,138)	139 (136,143)	<0.001
(60,70]	132 (129,136)	133 (129,136)	<0.001
(70,80]	134 (131,137)	134 (132,138)	<0.001
x15_20
(18,30]	105 (103,106)	105 (103,108)	<0.001
(30,40]	110 (108,112)	110 (108,112)	<0.001
(40,50]	111 (108,113)	111 (109,114)	<0.001
(50,60]	115 (112,118)	111 (108,115)	<0.001
(60,70]	114 (111,117)	109 (106,113)	<0.001
(70,80]	111 (108,115)	107 (103,110)	0.00141
x20_30
(18,30]	170 (166,174)	166 (162,172)	<0.001
(30,40]	182 (178,187)	172 (167,175)	<0.001
(40,50]	182 (177,187)	171 (166,177)	<0.001
(50,60]	186 (179,192)	168 (161,174)	<0.001
(60,70]	176 (171,180)	151 (147,158)	<0.001
(70,80]	136 (130,143)	118 (112,123)	0.00141
x30_40
(18,30]	95 (93,98)	87 (85,92)	0.02278
(30,40]	96 (92,99)	81 (78,85)	<0.001
(40,50]	89 (85,93)	75 (71,79)	0.00277
(50,60]	82 (77,86)	64 (59,68)	0.11019
(60,70]	67 (62,71)	45 (40,48)	0.79112
(70,80]	35 (32,38)	24 (22,27)	0.83662
x40_50
(18,30]	29 (28,31)	26 (24,27)	0.02278
(30,40]	28 (26,30)	22 (19,23)	<0.001
(40,50]	24 (22,26)	19 (17,21)	0.00277
(50,60]	20 (18,21)	13 (11,14)	0.11019
(60,70]	14 (12,16)	7 (6,8)	0.79112
(70,80]	5 (4,6)	3 (3,4)	0.83662
x50_60
(18,30]	6 (6,7)	6 (6,7)	0.89984
(30,40]	6 (6,7)	5 (4,6)	0.84073
(40,50]	5 (5,6)	4 (4,5)	0.75824
(50,60]	4 (3,4)	3 (2,3)	0.14261
(60,70]	2 (2,3)	2 (1,2)	0.03213
(70,80]	1 (0,1)	1 (0,1)	0.05500
x60_inf
(18,30]	3 (3,4)	4 (4,5)	0.89984
(30,40]	3 (2,3)	3 (3,4)	0.84073
(40,50]	3 (2,3)	3 (3,3)	0.75824
(50,60]	1 (1,2)	2 (2,3)	0.14261
(60,70]	1 (1,1)	1 (1,1)	0.03213
(70,80]	0 (0,0)	0 (0,0)	0.05500
total
(18,30]	14641 (14419,14863)	14423 (14202,14824)	0.21164
(30,40]	15034 (14712,15357)	14272 (13949,14528)	<0.001
(40,50]	14568 (14216,14920)	13865 (13513,14219)	0.00306
(50,60]	14102 (13758,14446)	13017 (12673,13341)	<0.001
(60,70]	12967 (12665,13269)	11428 (11127,11722)	<0.001
(70,80]	10433 (10131,10736)	9607 (9305,9873)	<0.001
mean
(18,30]	10 (10,10)	10 (10,10)	0.21164
(30,40]	10 (10,11)	10 (10,10)	<0.001
(40,50]	10 (10,10)	10 (9,10)	0.00306
(50,60]	10 (10,10)	9 (9,9)	<0.001
(60,70]	9 (9,9)	8 (8,8)	<0.001
(70,80]	7 (7,7)	7 (6,7)	<0.001

Code

pa_mims %>%
  left_join(demo_gh) %>% 
  filter(age_in_years_at_screening >= 18) %>% 
  pivot_longer(cols = c(total, mean))  %>%
  mutate(name = factor(name, labels = c("Mean MIMS", "Total MIMS"))) %>% 
  ggplot(aes(x = age_in_years_at_screening, y = value, color = gender)) +
  geom_jitter(alpha = .1, size = .7, width = .25) +
  geom_smooth(se = F) +
  facet_wrap(name~., scales = "free_y") + 
  labs(x = "Age", y = "Value", color = "")+
  scale_color_manual(values = c(col1, col2)) +
  theme(legend.position = "bottom",
        legend.title = element_blank())

Code

pa_mims %>%
  left_join(demo_gh) %>% 
  filter(age_in_years_at_screening >= 18) %>% 
  pivot_longer(cols = starts_with("x")) %>% 
  mutate(name = factor(name, levels = c("x0_2", "x2_5", "x5_10", "x10_15", "x15_20", "x20_30", "x30_40", "x40_50", "x50_60", "x60_inf"))) %>%
  ggplot(aes(x = age_in_years_at_screening, y = value, color = gender)) +
  geom_jitter(alpha = .1, size = .7, width = .25) +
  geom_smooth(se = F) +
  facet_wrap(name~., scales = "free_y") + 
  labs(x = "Age", y = "Minutes", color = "")+
  scale_color_manual(values = c(col1, col2)) +
  theme(legend.position = "bottom",
        legend.title = element_blank())

Code

survey_design =
  df_small_gh_mims %>%
  filter(age_in_years_at_screening >= 18) %>%
  mutate(weight = full_sample_2_year_mec_exam_weight / 2,
         weight_norm = weight / mean(weight)) %>%
  ungroup()  %>%
  survey::svydesign(
    id = ~ masked_variance_pseudo_psu,
    strata = ~ masked_variance_pseudo_stratum,
    weights = ~ weight_norm,
    data = .,
    nest = TRUE
  )


calc_by_age =
  function(var, df) {
    # var = "total_oaksteps"
    formula = as.formula(paste("~", var))
    total_by_age_gender = svyby(formula,
                                ~ age_in_years_at_screening + gender,
                                df,
                                svymean) %>%
      rename(mean = contains(var)) %>%
      mutate(metric = var)
  }


means_df =
  map_dfr(.x = pa_mims %>% select(starts_with("x"), total, mean) %>% colnames(),
          .f = calc_by_age, df = survey_design)

# get smooths for means and confidence intervals
models = means_df %>%
  mutate(lb = mean - 1.96 * se,
         ub = mean + 1.96 * se) %>%
  tidyr::nest(data = -c(metric, gender)) %>%
  dplyr::mutate(
    # Perform loess calculation on each group
    m = purrr::map(data, loess,
                   formula = mean ~ age_in_years_at_screening, span = .75),
    # Retrieve the fitted values from each model
    fitted_mean = purrr::map(m, `[[`, "fitted"),
    l = purrr::map(data, loess,
                   formula = lb ~ age_in_years_at_screening, span = .75),
    # Retrieve the fitted values from each model
    fitted_lb = purrr::map(l, `[[`, "fitted"),
    u = purrr::map(data, loess,
                   formula = ub ~ age_in_years_at_screening, span = .75),
    # Retrieve the fitted values from each model
    fitted_ub = purrr::map(u, `[[`, "fitted")
  )

# Apply fitted y's as a new column
results = models %>%
  dplyr::select(-m, -l, -u) %>%
  tidyr::unnest(cols = c(data, contains("fitted")))

results %>% 
  filter(!grepl("x", metric)) %>% 
  mutate(metric = factor(metric, labels = c("Mean MIMS", "Total MIMS"))) %>% 
    ggplot(aes(x = age_in_years_at_screening, y = fitted_mean,
             ymin = fitted_lb, ymax = fitted_ub, color = gender, fill = gender)) +
  facet_wrap(. ~ metric, scales = "free_y") +
  geom_line(linewidth = 1) +
  geom_ribbon(alpha = .2, linetype = 0) +
  theme_light() +
  theme(legend.position = "bottom",
        legend.title = element_blank(),
        text =element_text(size = 14),
        strip.text = element_text(size = 12),
        panel.grid.minor = element_blank()) +
  labs(x = "Age (years)", y = "Mean Value",
       title = "Smoothed survey weighted mean daily minutes by age")+
  scale_x_continuous(breaks=seq(20,80,10)) +
  scale_color_manual(values = c(col1, col2))+
  scale_fill_manual(values = c(col1, col2))

Code

results %>% 
  filter(grepl("x", metric)) %>% 
  mutate(metric = factor(metric, levels = c("x0_2", "x2_5", "x5_10", "x10_15", "x15_20", "x20_30", "x30_40", "x40_50", "x50_60", "x60_inf"))) %>% 
    ggplot(aes(x = age_in_years_at_screening, y = fitted_mean,
             ymin = fitted_lb, ymax = fitted_ub, color = gender, fill = gender)) +
  facet_wrap(. ~ metric, scales = "free_y") +
  geom_line(linewidth = 1) +
  geom_ribbon(alpha = .2, linetype = 0) +
  theme_light() +
  theme(legend.position = "bottom",
        legend.title = element_blank(),
        text =element_text(size = 14),
        strip.text = element_text(size = 12),
        panel.grid.minor = element_blank()) +
  labs(x = "Age (years)", y = "Mean Daily Minutes",
       title = "Smoothed survey weighted mean daily minutes by age")+
  scale_x_continuous(breaks=seq(20,80,10)) +
  scale_color_manual(values = c(col1, col2))+
  scale_fill_manual(values = c(col1, col2))

NYFFS

Tables

Code

demo_y = readRDS(here::here("data", "demographics", "subset_Y_tidy.rds")) 
pa_y = readRDS(here::here("data", "accelerometry", "summarized", 
                          "pa_summary_AC_paxy.rds")) %>% 
  mutate(SEQN = as.character(SEQN)) %>% 
  group_by(SEQN) %>%
  summarize(across(-PAXDAYM, ~mean(.x, na.rm = TRUE)))

df_small_y = pa_y %>% 
  left_join(demo_y) %>% 
  select(SEQN,
         gender,
         masked_variance_pseudo_psu,
         masked_variance_pseudo_stratum,
         full_sample_interview_weight,
         full_sample_mec_exam_weight,
         age_in_years_at_screening,
         total, 
         mean,
         contains("montoye"),
         starts_with("x")) 
labs = c("SEQN", "Gender",  
         "Pseudo PSU",
         "Psueudo stratum",
         "2 yr int weight", "2 yr exam weight",
         "Age", "Total AC", "Mean AC", "Sedentary minutes", "Light minutes",  "Moderate minutes", "Vigorous minutes",
         "MVPA minutes",
         "MVPA (10 min bout)",
         "Vig (10 min bout)",
         "Time <100", "Time [100,500]","Time (500,1000]","Time (1000,2000]", "Time (2000,3000]", "Time (3000,4000]", "Time (4000,5000]", "Time (5000,6000]", "Time (6000,7000]",  "Time (7000,8000]", "Time (8000,9000]", "Time >9000") 

names(labs) = colnames(df_small_y)

df_small_y =
  df_small_y %>%
  labelled::set_variable_labels(!!!labs)

df_svy =
  df_small_y %>% 
  select(SEQN,
         gender,
         masked_variance_pseudo_psu,
         masked_variance_pseudo_stratum,
         full_sample_interview_weight,
         full_sample_mec_exam_weight,
         age_in_years_at_screening,
         total, 
         mean,
         mvpa_montoye_bout,
         vig_montoye_bout,
         contains("range")) %>% 
  mutate(WTMEC4YR_norm = full_sample_mec_exam_weight/mean(full_sample_mec_exam_weight, na.rm = TRUE)) %>%
  svydesign(ids = ~masked_variance_pseudo_psu, weights = ~WTMEC4YR_norm,
            strata = ~masked_variance_pseudo_stratum, data=., nest=TRUE)

# survey weighted table
df_svy %>%
  tbl_svysummary(
    by = gender,
    include = -c(SEQN, masked_variance_pseudo_psu, masked_variance_pseudo_stratum, full_sample_interview_weight, full_sample_mec_exam_weight, WTMEC4YR_norm),
    statistic = list(
      all_continuous() ~ "{mean} ({sd})",
      all_categorical() ~ "{n} ({p}%)"
    ),
    digits = all_continuous() ~ 0,
    missing_text = "Missing",
  ) %>%
  add_overall() %>%
  add_p()

Characteristic	Overall N = 1,370 1	Female N = 674 1	Male N = 696 1	p-value 2
Age	9 (4)	9 (4)	9 (4)	0.2
Total AC	3,538,697 (805,142)	3,546,168 (755,892)	3,531,471 (850,561)	0.7
Mean AC	2,460 (559)	2,465 (525)	2,455 (591)	0.7
MVPA (10 min bout)	181 (85)	184 (83)	177 (88)	0.10
Vig (10 min bout)	80 (49)	79 (46)	81 (52)	0.7
1 Mean (SD)
2 Design-based KruskalWallis test

Code

survey_design =
  df_small_y %>%
  mutate(age_cat = cut(age_in_years_at_screening,
                       breaks=c(3, 6, 9, 12, 15))) %>%
  mutate(weight = full_sample_mec_exam_weight,
         weight_norm = weight / mean(weight)) %>%
  ungroup()  %>%
  survey::svydesign(
    id = ~ masked_variance_pseudo_psu,
    strata = ~ masked_variance_pseudo_stratum,
    weights = ~ weight_norm,
    data = .,
    nest = TRUE
  )

calc_by_age =
  function(var, df) {
    formula = as.formula(paste("~", var))
    total_by_age_gender = svyby(formula,
                                ~ age_cat + gender,
                                df,
                                svymean) %>%
      rename(mean = contains(var)) %>%
      mutate(metric = var)
  }

means_df =
  map_dfr(.x = pa_y %>% select(contains("bout"), contains("x"), contains("montoye")) %>% colnames(),
          .f = calc_by_age, df = survey_design)


small = 
  means_df %>%
  pivot_wider(names_from = gender, values_from = c(mean, se)) %>% 
  mutate(val_female = paste0(round(mean_Female), " (", round(mean_Female - 1.96 * se_Female), ",", round(mean_Female + 1.96 * se_Female), ")"),
         val_male = paste0(round(mean_Male), " (", round(mean_Male - 1.96 * se_Female), ",", round(mean_Male + 1.96 * se_Male), ")"))

# test = svyttest(mvpa_montoye ~ gender, survey_design)
pvals = c()
for(i in 1:nrow(small)){
  row = means_df[i,]
  cat = row$age_cat
  var = row$metric
  survey_design =
    df_small_y %>%
    mutate(age_cat = cut(age_in_years_at_screening,
                       breaks=c(3, 6, 9, 12, 15))) %>%
    filter(age_cat == cat) %>%
    mutate(weight = full_sample_mec_exam_weight,
           weight_norm = weight / mean(weight)) %>%
    ungroup()  %>%
    survey::svydesign(
      id = ~ masked_variance_pseudo_psu,
      strata = ~ masked_variance_pseudo_stratum,
      weights = ~ weight_norm,
      data = .,
      nest = TRUE
    )
  formula = as.formula(paste(var,  "~ gender"))
  test = svyttest(formula, survey_design)
  pvals = append(pvals, unname(test$p.value))
}

small %>% 
  select(metric, age_cat, Female = val_female, Male = val_male, mean_Female, mean_Male) %>% 
  bind_cols(p = pvals) %>% 
  mutate(pval = if_else(p < 0.001, "<0.001", format.pval(p, digits = 2)))  %>%
  group_by(metric) %>% 
  gt() %>% 
  cols_hide(c("mean_Female", "mean_Male", "p")) %>% 
  tab_style(
    # Style to apply when conditions are met
    style = cell_fill(color = "lightpink"),
    locations = cells_body(
      rows = mean_Female > mean_Male & p < 0.05
    )
  ) %>% 
  tab_style(
    # Style to apply when conditions are met
    style = cell_fill(color = "lightblue"),
    locations = cells_body(
      rows = mean_Female < mean_Male & p < 0.05
    )
  )

age_cat	Female	Male	pval
mvpa_montoye_bout
(3,6]	243 (227,258)	233 (218,252)	0.16257
(6,9]	213 (193,232)	193 (174,205)	0.03998
(9,12]	159 (148,170)	153 (141,162)	0.30361
(12,15]	115 (100,130)	108 (93,118)	0.41696
vig_montoye_bout
(3,6]	103 (93,113)	107 (97,117)	0.16257
(6,9]	96 (83,109)	92 (80,102)	0.03998
(9,12]	68 (61,74)	70 (64,75)	0.30361
(12,15]	47 (39,54)	49 (42,55)	0.41696
x0_100
(3,6]	563 (555,571)	551 (543,560)	0.42823
(6,9]	543 (533,553)	533 (523,540)	0.58236
(9,12]	551 (544,559)	539 (531,548)	0.49074
(12,15]	578 (557,599)	581 (561,600)	0.55775
x100_500
(3,6]	67 (64,69)	72 (70,76)	0.42823
(6,9]	73 (69,77)	84 (80,88)	0.58236
(9,12]	91 (85,97)	106 (101,112)	0.49074
(12,15]	114 (109,118)	130 (125,135)	0.55775
x500_1e_03
(3,6]	65 (63,68)	72 (70,75)	0.02874
(6,9]	76 (73,79)	86 (82,89)	0.13839
(9,12]	93 (88,97)	104 (99,107)	0.06680
(12,15]	110 (107,113)	121 (118,125)	0.76776
x1e_03_2e_03
(3,6]	113 (109,118)	123 (118,127)	0.02874
(6,9]	132 (126,138)	142 (136,147)	0.13839
(9,12]	152 (146,158)	157 (152,160)	0.06680
(12,15]	160 (156,164)	165 (161,170)	0.76776
x2e_03_3e_03
(3,6]	108 (103,112)	110 (105,115)	0.01209
(6,9]	116 (110,121)	119 (113,122)	0.00154
(9,12]	118 (115,121)	116 (113,119)	<0.001
(12,15]	111 (108,115)	106 (103,110)	<0.001
x3e_03_4e_03
(3,6]	102 (98,105)	102 (98,105)	0.01209
(6,9]	101 (97,106)	101 (96,104)	0.00154
(9,12]	96 (94,97)	92 (90,95)	<0.001
(12,15]	86 (83,89)	79 (76,82)	<0.001
x4e_03_5e_03
(3,6]	94 (91,96)	91 (89,93)	0.00191
(6,9]	88 (85,91)	85 (82,88)	0.00159
(9,12]	80 (77,83)	76 (73,78)	0.00121
(12,15]	70 (67,73)	63 (61,66)	<0.001
x5e_03_6e_03
(3,6]	81 (79,83)	78 (75,80)	0.00191
(6,9]	75 (73,77)	69 (67,71)	0.00159
(9,12]	66 (63,68)	62 (60,65)	0.00121
(12,15]	58 (54,61)	52 (49,55)	<0.001
x6e_03_7e_03
(3,6]	66 (65,68)	63 (61,65)	<0.001
(6,9]	61 (58,63)	55 (52,57)	0.01115
(9,12]	53 (51,56)	50 (47,52)	0.07327
(12,15]	46 (43,49)	42 (39,44)	0.06065
x7e_03_8e_03
(3,6]	51 (49,53)	47 (45,50)	<0.001
(6,9]	47 (44,49)	42 (39,43)	0.01115
(9,12]	41 (39,43)	37 (35,39)	0.07327
(12,15]	34 (32,37)	30 (27,32)	0.06065
x8e_03_9e_03
(3,6]	37 (35,39)	35 (33,37)	0.27116
(6,9]	35 (32,37)	30 (28,31)	0.19712
(9,12]	30 (29,32)	27 (25,29)	0.44564
(12,15]	24 (22,27)	20 (18,22)	0.04788
x9e_03_inf
(3,6]	91 (84,99)	97 (90,104)	0.27116
(6,9]	93 (84,101)	94 (85,102)	0.19712
(9,12]	68 (62,74)	71 (65,76)	0.44564
(12,15]	47 (41,53)	48 (41,53)	0.04788
sed_montoye
(3,6]	901 (890,913)	913 (901,928)	0.95553
(6,9]	925 (909,941)	948 (932,958)	0.70090
(9,12]	990 (973,1006)	1008 (992,1024)	0.03792
(12,15]	1059 (1037,1082)	1091 (1069,1109)	0.00855
light_montoye
(3,6]	111 (107,115)	111 (107,115)	0.95553
(6,9]	111 (106,116)	111 (106,115)	0.70090
(9,12]	105 (103,108)	102 (100,105)	0.03792
(12,15]	95 (92,98)	87 (84,91)	0.00855
mod_montoye
(3,6]	151 (147,155)	146 (142,149)	0.10790
(6,9]	141 (137,145)	134 (130,139)	0.11076
(9,12]	127 (122,131)	120 (116,125)	0.02906
(12,15]	111 (106,116)	101 (96,105)	0.00849
vig_montoye
(3,6]	276 (264,288)	270 (258,284)	0.10790
(6,9]	262 (246,278)	246 (230,255)	0.11076
(9,12]	217 (205,229)	208 (196,217)	0.02906
(12,15]	172 (157,187)	159 (144,170)	0.00849
mvpa_montoye
(3,6]	472 (459,484)	461 (448,476)	0.04812
(6,9]	447 (430,463)	424 (408,434)	<0.001
(9,12]	385 (369,401)	368 (352,382)	0.03768
(12,15]	320 (300,341)	293 (273,309)	0.04173

Smoothed, unweighted

Code

pa_y %>%
  left_join(demo_y) %>% 
  pivot_longer(cols = contains("bout")) %>% 
  mutate(name = factor(name, labels = c("MVPA", "VPA"))) %>% 
  ggplot(aes(x = age_in_years_at_screening, y = value, color = gender)) +
  geom_jitter(alpha = .1, width = .25) +
  geom_smooth(se = F) +
  facet_grid(name~., scales = "free_y") + 
  labs(x = "Age", y = "Minutes in MVPA", color = "")+
  scale_y_continuous(breaks=seq(0,600,60)) + 
  scale_color_manual(values = c(col1, col2)) + 
  theme(legend.position = c(0.9, 0.9),
        legend.title = element_blank()) + 
  scale_x_continuous(breaks=seq(3,15,2))

Smoothed, weighted

Code

survey_design =
  df_small_y %>%
  mutate(weight = full_sample_mec_exam_weight,
         weight_norm = weight / mean(weight)) %>%
  ungroup()  %>%
  survey::svydesign(
    id = ~ masked_variance_pseudo_psu,
    strata = ~ masked_variance_pseudo_stratum,
    weights = ~ weight_norm,
    data = .,
    nest = TRUE
  )


calc_by_age =
  function(var, df) {
    # var = "total_oaksteps"
    formula = as.formula(paste("~", var))
    total_by_age_gender = svyby(formula,
                                ~ age_in_years_at_screening + gender,
                                df,
                                svymean) %>%
      rename(mean = contains(var)) %>%
      mutate(metric = var)
  }

means_df =
  map_dfr(.x = pa_y %>% select(contains("bout")) %>% colnames(),
          .f = calc_by_age, df = survey_design)

# get smooths for means and confidence intervals
models = means_df %>%
  mutate(lb = mean - 1.96 * se,
         ub = mean + 1.96 * se) %>%
  tidyr::nest(data = -c(metric, gender)) %>%
  dplyr::mutate(
    # Perform loess calculation on each group
    m = purrr::map(data, loess,
                   formula = mean ~ age_in_years_at_screening, span = .75),
    # Retrieve the fitted values from each model
    fitted_mean = purrr::map(m, `[[`, "fitted"),
    l = purrr::map(data, loess,
                   formula = lb ~ age_in_years_at_screening, span = .75),
    # Retrieve the fitted values from each model
    fitted_lb = purrr::map(l, `[[`, "fitted"),
    u = purrr::map(data, loess,
                   formula = ub ~ age_in_years_at_screening, span = .75),
    # Retrieve the fitted values from each model
    fitted_ub = purrr::map(u, `[[`, "fitted")
  )

# Apply fitted y's as a new column
results = models %>%
  dplyr::select(-m, -l, -u) %>%
  tidyr::unnest(cols = c(data, contains("fitted")))

results %>%
  mutate(metric = factor(metric, labels = c("MVPA", "VPA"))) %>% 
    ggplot(aes(x = age_in_years_at_screening, y = fitted_mean,
             ymin = fitted_lb, ymax = fitted_ub, color = gender, fill = gender)) +
  facet_grid(. ~ metric) +
  geom_line(linewidth = 1) +
  geom_ribbon(alpha = .2, linetype = 0) +
  theme_light() +
  theme(legend.position = "bottom",
        legend.title = element_blank(),
        text =element_text(size = 14),
        strip.text = element_text(size = 12),
        panel.grid.minor = element_blank()) +
  labs(x = "Age (years)", y = "Mean Daily Minutes",
       title = "Smoothed survey weighted mean daily minutes by age ")+
  scale_x_continuous(breaks=seq(3, 15, 2)) +
  scale_y_continuous(breaks=seq(0,600,30)) + 
  scale_color_manual(values = c(col1, col2))+
  scale_fill_manual(values = c(col1, col2))

Code

means_df =
  map_dfr(.x = pa_y %>% select(starts_with("x")) %>% colnames(),
          .f = calc_by_age, df = survey_design)

# get smooths for means and confidence intervals
models = means_df %>%
  mutate(lb = mean - 1.96 * se,
         ub = mean + 1.96 * se) %>%
  tidyr::nest(data = -c(metric, gender)) %>%
  dplyr::mutate(
    # Perform loess calculation on each group
    m = purrr::map(data, loess,
                   formula = mean ~ age_in_years_at_screening, span = .75),
    # Retrieve the fitted values from each model
    fitted_mean = purrr::map(m, `[[`, "fitted"),
    l = purrr::map(data, loess,
                   formula = lb ~ age_in_years_at_screening, span = .75),
    # Retrieve the fitted values from each model
    fitted_lb = purrr::map(l, `[[`, "fitted"),
    u = purrr::map(data, loess,
                   formula = ub ~ age_in_years_at_screening, span = .75),
    # Retrieve the fitted values from each model
    fitted_ub = purrr::map(u, `[[`, "fitted")
  )

# Apply fitted y's as a new column
results = models %>%
  dplyr::select(-m, -l, -u) %>%
  tidyr::unnest(cols = c(data, contains("fitted")))

results %>% 
    ggplot(aes(x = age_in_years_at_screening, y = fitted_mean,
             ymin = fitted_lb, ymax = fitted_ub, color = gender, fill = gender)) +
  facet_wrap(. ~ metric, scales = "free_y") +
  geom_line(linewidth = 1) +
  geom_ribbon(alpha = .2, linetype = 0) +
  theme_light() +
  theme(legend.position = "bottom",
        legend.title = element_blank(),
        text =element_text(size = 14),
        strip.text = element_text(size = 12),
        panel.grid.minor = element_blank()) +
  labs(x = "Age (years)", y = "Mean Daily Minutes",
       title = "Smoothed survey weighted mean daily minutes by age ")+
    scale_x_continuous(breaks=seq(3, 15, 2)) +
  scale_color_manual(values = c(col1, col2))+
  scale_fill_manual(values = c(col1, col2))